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Preface 


Many problems in mathematical physics and engineering are modeled by elliptic 
systems of partial differential equations. Well-known examples in this context are 
higher-dimensional steady-state heat conduction [11], acoustics (see [8, 15]), grav- 
itational potential [10], fluid mechanics [14], and elasticity theory (plane strain, 
bending of thin plates, and stationary flexural oscillations—see [11, 12]). Among 
the solution techniques for such problems, a prominent role is played by boundary 
integral equation methods (BIEMs), which, apart from being powerful and elegant, 
have some decisive advantages over other procedures, an important one being that 
they change a problem from its formulation in terms of an unbounded differential 
operator to one for an integral operator, thus making it more appealing and tractable 
from an analytic viewpoint. In essence, the ellipticity of the problems is shifted to 
the boundary, where it gives rise to integral equations that are then solved in suit- 
able function spaces [6]. BIEMs also yield closed form solutions, which renders 
them very useful for numerical computation. 

The many different types of BIEMs developed so far [6] can be divided into 
two main categories: direct and indirect. The former employ what might be called 
‘designer’ solutions, of a form chosen purely for mathematical convenience. By 
contrast, the unknown functions in the latter have physical significance; in elasticity, 
for example, they may characterize the displacement field or the force vector acting 
on the boundary of the domain. 

An elastic body under the action of external forces experiences deformation, 
mathematically described by means of the stress and strain tensors and the dis- 
placement vector (see [9, 1]). These quantities satisfy a system of equations that, 
under certain conditions, can be reduced to one in only two—rather than the orig- 
inal three—independent space variables. Such is the case of the system governing 
the equilibrium of thin elastic plates with transverse shear deformation, extensively 
investigated in [6, 7] in spaces of smooth functions, and in [2] in spaces of distribu- 
tions. The corresponding stationary oscillations model has been studied in [16], and 
that of dynamic deformations in [3]. 

In this book, we consider the system of equations (known as the Winkler model 
[18]) that describes the equilibrium of a thin elastic plate with in-plane deformation 
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and no bending, which lies on an elastic foundation and is subjected to Dirichlet, 
Neumann, or Robin boundary conditions. This model has many important applica- 
tions in engineering problems arising in geotechnical research, road construction, 
biomechanics, and other practical fields. A brief presentation of some preliminary 
results can be found in [4, 5]. Our intention is to describe the mathematical model 
analytically and then use it to show how a boundary element method, based on 
the boundary integral equation technique, can be constructed and manipulated to 
compute an approximate (numerical) solution. The advantage of this type of ap- 
proach over the use of finite elements or other classical computational procedures is 
twofold: it reduces the original two-dimensional setup to a one-dimensional prob- 
lem, and provides a faster rate of convergence. 

The material in the book is organized in five chapters. 

In Chapter 1, we describe the model and list the main boundary value problems 
associated with its governing system. 

Chapter 2 contains the definition of the elastic potentials in terms of a matrix of 
fundamental solutions, and a description of the mapping properties of the boundary 
integral operators generated by these potentials. 

The main thrust in Chapter 3 is the use of the layer potentials in the solution of our 
boundary value problems by means of direct and classical indirect BIEMs. These 
results are fundamental in establishing the well-posed nature of the mathematical 
model. 

In Chapter 4, we make a detailed presentation of the Mathematica® software in 
relation to its use in our numerical handling of the problems. 

Finally, Chapter 5 consists of a collection of computational examples that illus- 
trate the implementation of the collocation method for the direct and classical indi- 
rect BIEMs with various conditions prescribed on smooth and non-smooth bound- 
aries, and different choices of splines. 

The book should be a good source of information for readers who want to get 
at-a-glance theoretical and practical details about the analytic structure and numeri- 
cal applications (in the Mathematica® environment) of boundary integral equation 
methods. 


We wish to express our thanks to Elizabeth Loew, Executive Editor for Math- 
ematics at Springer, who has guided the evolution of this project with impeccable 
professionalism and great efficiency. 

Also, two of us (CC and DD) are grateful to our wives for their support, under- 
standing, and remarkable staying power in the face of stiff competition from our 
computers during the production of this book. 
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Chapter 1 
The Mathematical Model 


1.1 Basic Equations 


Throughout the book, Latin and Greek subscripts take the values 1,2,3 and 1,2, 
respectively, and the convention of summation over repeated indices is understood. 
For simplicity, we denote by J both the identity matrix on any space of square ma- 
trices and the identity operator on any space of functions. Also, we denote the trans- 
pose of a matrix X by X" and the derivatives of a function f = f(x;) by 
of 
—_ = i 0; ; 
Ox; f ? if 
af 
a = fip= dif, 
OxjOXx j f J if 


with the obvious generalization for higher-order derivatives. 
Let S be a finite domain in R?, with a simple and closed boundary 0S oriented in 
the positive direction. 


1.1 Definition. (i) A function f defined in S (on OS) is called Hélder continuous 
with index a € (0,1) if there is c = const > 0 such that 


If(x) -—fO)| <elx—yl™ Vx,y ES (x,y € OS). 


The space of all such functions is denoted by C°-™(S) (C°(dS)). 


(ii) A function f defined in S (on OS) is called Hélder continuously differen- 
tiable with index a € (0,1) if it is differentiable in S (on 0S) and its derivatives 
belong to C°(S) (C°“(AS)). The space of all such functions is denoted by C!:“(S) 
(c!(dS)). 


1.2 Remark. If S is an infinite domain, then a function f defined on S is Hélder 
continuous (Hélder continuously differentiable) in S if it has this property on any 
finite subdomain of S. Of course, the constant c in Definition 1.1 will vary with the 
subdomain. 
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2: 1 The Mathematical Model 


We now assume that the boundary curve 0S is of class C?, denote S by S+, and 


write 
S- =R*\(StUas), St=S*tUaS, 5 =S- VAS. 


Let x = (x1,x2) be a generic point in a Cartesian system of coordinates with the 
origin (0,0) € St, and let v(x) be the outward unit normal at x on OS (see Fig. 1.1). 
For x, y € R*, we write 


Ix —yl = [1 —y1)? + (2-27. 


x2 


v(x) 


os 


x] 


Fig. 1.1 The geometric configuration of the domains and boundary. 


The equilibrium equations of three-dimensional elasticity are 
tijg+qi =9, 


where the f;; = tj; are the components of the stress tensor and q; are the components 
of the body force vector. For a homogeneous and isotropic material, the constitutive 
equations have the form 


tij= Aug .0ij + M(ui,j + Uji), 


where u; are the displacement components, A and wu are the elastic (Lamé) coef- 
ficients, and 6; ; are the Kronecker delta. The components of the stress vector in a 
direction n = (nj,n2,n3)" are 

B= jj, 


and the internal energy density is 


1.1 Basic Equations 3 


The model of plain strain is based on the assumption that 


n= (n,n2,0)'. 


Hence, the equilibrium equations, the constitutive relations, the stress tensor com- 
ponents, and the internal energy density can be expressed as 


tap.p + qa = 0, (1.1) 
top = Auy,yOop + U(Ua,p + Upc)» (1.2) 

ta =tapng; (1.3) 
Baits. (1.4) 


Combining (1.1) and (1.2), we find that 


(A+ )uB Ba + HUa,pp +a = 9, 
or 
(A + W)grad divu+ pAu = q, (1.5) 
which can be written as 
A(O1,02)u=4q, 
where 
u= (u,u2)', q=(41,92)', A = 07+, 


and A(d1, 02) is the matrix differential operator 


atana= (Mes (A +p) aa ) 


(A+M)A0, bA+(A+H) dz 
In the homogeneous case—that is, in the absence of body forces—system (1.5) 


becomes 
A(a,,)u=0. 


If the plate lies on an elastic foundation, the differential operator changes to 


oy. 
uA+(Atp)ae—k = (A+p)a1a i 4a 


Z(A1,) = 
(1,02) ( (A+M)A® pwA+(A+p)a2—-k 


or 
Z(O1,02) =A(A, 02) — KI, (1.7) 


where k = const > 0 is the elastic constant of the foundation material. Then the 
homogeneous equilibrium system can be written as 


4 1 The Mathematical Model 
Z(0,02)u=0. (1.8) 


1.3 Theorem. /f 
A+pu>0, u>O, (1.9) 


then system (1.8) is elliptic. 


Proof. Since, under conditions (1.9), the system A(Q),02)u = 0 is elliptic (see The- 
orem 2.1 in [6]), from (1.7) it follows that so is (1.8). 


1.4 Remark. In what follows, we assume that inequalities (1.9) are satisfied. 


Combining (1.2) and (1.3) and setting ng = Vq (the components of the outward 
unit normal to 0S), we arrive at 


to, = top Vp = Tapup, 


where the boundary stress operator T (01, 02) is given by 


A+2 fa) O 0, +Avja 
T(a1,d%) = ( +2U)ViQ}+MV202 = L201 + AV 02 ) (1.10) 


Av20, + UV} 02 Lyd) +(A +2) v2.02 


From (1.2) and (1.4) it follows that the internal energy density in the absence of 
a foundation can be written as 


€=E(u,u) = 5 [(A + 2u) (uy + U5») +2Auj 142.2 + L(u12+ u2,1)’]. 
For the elastic foundation problem, the internal energy density is 
E=E(u,u) =@(u,u)+4ku'u. (1.11) 
1.5 Theorem. E(u, u) is a positive definite quadratic form. 


This assertion follows from Theorem 2.2 in [6], which states that &(u,u) is pos- 
itive definite, and the fact that the second term in (1.11) is also positive definite. 


1.2 Boundary Value Problems 


Let .@xq be the vector space of p x q matrices, and let ./ be the vector space of 
functions u € 4%, in S~ such that, in polar coordinates, 


u(r,@) =O(r-'-*), a>0, as re. (1.12) 


Let 7,2, 2, %, XH, L © C(AS) be given 2 x 1 matrix-valued functions, and 
let o € C(OAS) be a positive definite 2 x 2 matrix-valued function. We consider the 
following six basic boundary value problems: 


1.2 Boundary Value Problems 5 
(i) Interior Dirichlet problem (Dt): Find u € C?(S+)C!(S+) such that 
(Zu)(x) =0, xESt, u(x) = A(x), xe OS. 
(ii) Exterior Dirichlet problem (D~ ): Find u € C?(S~) NC!(S~) 9. such that 
(Zu)(x)=0, xES”, u(x) =A(x), xe OS. 
(iii) Interior Neumann problem (N+): Find u € C?(S+) NC!(S*) such that 
(Zu)(x) =0, xE St, (Tu)(x)= Qa), xE OS. 
(iv) Exterior Neumann problem (N~ ): Find u € C?(S~) NC!(S~)A.& such that 
(Zu)(x) =0, xES”, (Tu)(x) =A (x), xe OS. 


(v) Interior Robin problem (R*): Find u € C?(S+) VC! (S*) such that 


(Zu)(x) =0, xES™, (Tu+ou)(x)= 4 (x), xe OS. 
(vi) Exterior Robin Problem (R~): Find u € C?(S~) NC!(S~).@ such that 
(Zu)(x) =0, xES”, (Tu—ou)(x)= L(x), x€ OS. 


Any function that satisfies pointwise one of the above systems of equations is 
termed a regular solution, or, simply, a solution of that boundary value problem. 
The following assertion is analogous to Green’s first identity [13]. 


1.6 Theorem. [fu € C?(S+)MC!(S*), then 


[u'(2u)aa= ft (tuyds—2 [ E(u,u)da. 
st 


St os 


Proof. From (1.7) and Theorem 2.7 in [6] it follows that 


j[u'(2u)da= [@a-kyuda= fu (Au)aa— fu (hu) da 


St 


St St St 
= [ub(tu)ds—2 [ (uu)\da—k f uTuda 
os St SF 


= fut(rwyas—2( fstuyaa+ $e [Tuda) 
os St St 


= ft tuyas—2 f E(uu)da. 
os St 


The above assertion has an important consequence. 


6 1 The Mathematical Model 
1.7 Corollary (Betti formulas). (i) [fu € C?(S+) NC!(S*) and Zu = 0 in S*, then 
2 [ E(u,u)da= [ w"(Tu) as. (1.13) 
St os 
(ii) fu € C?(S~)NC!(S-) NL and Zu=0 in S~, then 
2 [ 8(u.u)da=— fu" (Tuas. (1.14) 
s- as 
Proof. (i) By Theorem 1.6, 
[(Gu)aa = (rw ds—2 f E(u,w) da, 
st as st 


and the result follows from the fact that Zu = 0 in S*. 


(ii) Let Kr be a disc centered at the origin and of radius R sufficiently large so 
that St C Kp (see Fig. 1.2). 


x2 


os, 


x} 


aa 


K,\S 


Fig. 1.2. The domain Kp \ S*. 


By (1.7), 


i u' (Zu) da = i: u' (Au) da—k / u'uda. 
Kr\S* Kr\S+ Kr\S+ 
Since Zu = 0 in Kr \ S* C S~, it follows that 


0= / u' (Au) da—k uluda. 


Kpr\S* Kpr\S* 
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Then, by Corollary 2.8 in [6], 


j= i ul (Tu) ds—2 7 &(u,u)da—k / u'uda. (1.15) 
o(Kr\S*) Kp\St Kp\St 
Since here the outward normal to 0S is directed into St (see Fig. 1.2), we have 
i u'(Tu)ds = / u' (Tu)ds — [e(twas, 
8(Kp\5+) aKp as 


which, when substituted in (1.15), yields 


0= [@rwds— fu"(ru)as 


OKr os 
= if &(u,u)da—k : ul uda. (1.16) 
Kp\S* Kr\S* 
By (1.12), 
u=O(r}*), a>O, 
so 
wu=w+u=O(R? 7%) in Kpyg+5 

therefore, 

uluda = O(R > 7%)O(R*) = O(R- 2%) +0 asR>&. (1.17) 

Kr\St 

Since 


Olu, Onu = Or), 
it follows that 


/ u" (Tu) ds = O(R7!-®)0(R-2-“)0(R) 


OKr 5 
=O(R-2-2%) +0 asR3o. (1.18) 


Also, Kr \ St expands into the whole of S- as R— ©; therefore, in the limit, for- 
mulas (1.16)—(1.18) generate (1.14). 


The analog of Green’s second identity [13] is called the reciprocity relation. 


1.8 Theorem. For any u, v € C?(S*)NC!(S*), 


/ [ul (zv) —v" (Zu) | da = / ul (Tv) —v"(Tu)]ds. 


St as 


8 1 The Mathematical Model 
Proof. By Theorem 2.9 in [6], 
/ tu? (Zv)—v"(Zu)] da = / u!(A—k)v—v" (A —Eda 


St St 


- u"(Av) — v¥(Au)]da~k [ (uTv— vu) da 


St St 
= 7 ul (Av) —v"(Aw)|da= / lu" (Tv) —v"(Tu)]|ds. 
St os 


1.9 Theorem (Uniqueness). Each of (D*), (D~), (N*), (N~), (R™), and (R7) has 
at most one solution. 


Proof. Let u= v—w, where v, w € C?(S+)MC!(S*) are any two solutions of one 
of (D*) or (N*) with the same nonhomogeneous term. Then u is a solution of the 
corresponding homogeneous problem, so, by (1.13), 


2 [ E(uu)da= | w"(Tu)as=0. 
st 


os 


Since E(u, u) is a positive definite quadratic form, it follows that E(u,u) = 0, which 
implies that u = 0, or v = w; hence, each of (D*) and (N*) has at most one solution. 

The proof for (D7) and (N ) is similar, with (1.14) used instead of (1.13). 

We now turn to the Robin problems. As above, the difference u = v— w of any 
two solutions v, w € C?(S*+)NC!(S*) of (R*) for a given .% and a given positive 
definite 2 x 2 matrix function o € C(0S), is a solution of the homogeneous problem 
(R*), so, again, from (1.13) it follows that 


2 | E(u,u)da= fe twas= ~ [ u(ou)ds, 
St 


as os 


or, what is the same, 


2 [ E(u,u)da+ | u™(ou)ds =. 
St os 


We have 
ulou= Ooplalps 
consequently, since E(u,u) is a positive definite quadratic form and o is a positive 
definite matrix, we conclude that u is a rigid displacement in S* that vanishes on 
OS, so u = 0, or v= w, which means that (Rt) has at most one solution. 
The proof for (R~) is similar to that for (R*), with the Betti formula (1.14) in- 
stead of (1.13). 


Chapter 2 
The Layer Potentials 


2.1 Fundamental Solutions 


We need to construct a matrix of fundamental solutions for the operator Z(0),02). 
This is a two-point 2 x 2 matrix function satisfying 


Z(0x)D(x,y) = —6(|x— yI)I, (2.1) 


where 0, indicates that differentiation is applied with respect to the point x and 6 is 
the Dirac delta distribution. If the scalar two-point function t(x,y) is a solution of 


the equation 
(detZ)(dr)t(x,y) = —6(|x—yI), 


then 
Z(Ox)|(adj Z)(Ax)t(x,y)I] = (detZ) (Ax )t(x, yl = —6(|x— yl), 
so 
D(x,y) = (adjZ)(x)t(x,y)I (2.2) 
is a matrix of fundamental solutions for Z. 
We have 
A+(A 0? —k A O10: 
(detz)(a,) = |4 (A+) dj (A+) z 
(A+M)O10, pUwA+(A+p)ds—k 
= p(A+2p)A* —2kpA —(A+p)kA +k? 
k(A +3) ke ) 
=p(A+2p)( A? 
+20) (a8 yd aoa 
= w(A+2u)(A —C7)(A -C5), (2.3) 
where ‘ 
Dr Date 
Ci > u ’ C5 a +2uU : 
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A fundamental solution for the operator A — hI, h = const, is [17] 
—(1/(22))Ko(hlx—yI), 


where Ko is the modified Bessel function of the second kind and order zero (see Fig. 
2.1); that is, 
(A —h*)Ko(h|x— y|) = —228(|x— yl). 


Ko(r) 


Fig. 2.1 The modified Bessel function Ko. 


This suggests that we seek t(x,y) of the form 
t(x,y) = pKo(Cilx—y|) +qKo(C2lx—yl), p,q = const, 


which then leads to 


= —2np6(|x—yl) + p(C{ —C3)Ko(Cilx—y|) — 26 (|x—y)). 
To eliminate 6, we take p = —q; hence, 
(A —C3)t(x,y) = p(Ct — Cz) Ko(Cilx—y)), 


so 
(A —Cf)(A — C3)t (x,y) = (A — CP) p(Ct — C2) Ko(Ci x — yl) 


= p(C} —C3)(A — CP) Ko(C\|x—y]) 
= —2np(Cj —Cz)5(|x—y)). (2.4) 
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Using (2.3) and (2.4), we find that 


(detZ) (0)¢(x,y) = M(A +2)(A —C})(A —C3)e(x,9) 
= —2mp(A +2 )A(C}-C3)4(|x—yI) 


k k 
= —2mpu(a+2n) (5 — 5 )athx—y) 
= —2mpk(A + u)6(|x—yI) 
=—6(|x—yI), 
where 
= 1 
P Inka +B)’ 
—* i . 
1 Fak +B)’ 
consequently, 
1 
Y) = Feta py KolCie— yl) — KolCalx—yI)}. (2.5) 


In this case, 
D(x,y) = (adjZ)(0,)t(x,y)I 
_ ("" + (A+p)02 —k (A + U)d;02 es 0 ) (2.6) 


—(A+p)da,  pA+(A+p)d? —k 0 t(x,y) 
Or, component-wise, 
Di (x,y) = HAt(x,y) Tv (A + 1) d3t(x,y) —kt(x,y), 
D12(x,y) = Dai (x,y) = —(A +H) oi ont (x,y), (2.7) 


Dro (x,y) = WAt(x,y) + (A + HOF (x,y) — (x,y), 


with all differentiation performed with respect to the point x. 
We notice the symmetry 


D(x,y) = D(y,x) =D" (x,y). (2.8) 


Another necessary tool for handling our boundary value problems is the matrix 
of singular solutions 
P(x,y) = [T(%)DV.2))", (2.9) 


where 
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Pi (X,¥) = H(A +2 )dyyyA (yt, y) +2U(A +H) Asyyt,12 29) 
—k[vi(y)(A + 2p )t,1 (x,y) + v2(y)t,2 Oy), 
Pia(x,y) = HOsyy A (y)t(x,y) — 2M (A + L)Os(yyty11 (x,y) 
—k[vi(y)t,2 xy) + va(y)t1 @y)], 
P21 (X,Y) = ~My A(y)t(x,y) +2U(A + W)A(yyt,22 Y) 
Vi (Y)t,2 (4,9) + Voy)t1 ey), 
Poo(X,¥) = MA + 2M) dy A (yt, y) — 2U(A +H) Asyyt,12 9) 
— [vi (y) (A + 2p )t,1 (x,y) — V20)t,2 Oy); 


here, we have denoted the normal and tangential derivatives at y on 0S by 


(2.10) 


= 


2 
) 


yy) = 9/AV(y) = vily) a1 + V2(y)o2, 
O.(y) = 9 /Os(y) = —V2(y)O1 + vi (y)o2, 
with differentiation in terms of the point y. 


2.1 Theorem. For x # y, the columns D‘“ of D and the columns P of P are 
solutions of (1.8). 


Proof. The first part of the assertion follows from the definition of D in (2.1). In 
fact, we can write 
Z(Ox)D(x,y) =0, x#y. (2.11) 


Furthermore, in view of the symmetry (2.8), we also have 
Z(dy)D(x,y) = 0, Z(0,)D(y,x) = 0, x#Y. 


Finally, by (2.11) and (2.9), we see that for x £ y, 


It is important to know the asymptotic behavior of D(x,y) and P(x,y) for small 
and large values of |x — y|. 
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First, for x close to y we have [6] 
Ko(§)=—(1+ 467? + Sit) Ing, 
so 


Ko(Ci|x — y|) — Ko(C2|x — y]) 


= 4(C-Ci)r in|x—yl + (CG -Ch)be—yl* nx y| + 
+4 (Con, — yy InC))|x—y|? + | (CG In —CflnC;)|x—y|* +---, 
and, by (2.5), 
t(x,y) = Alx—yl?In(|x— yl) + Blx—y]? + O(|x—yl*In x — yl), (2.12) 


where A and B are combinations of A and . Direct computation now shows that, in 
terms of the point x, 


Poh 
|x—y| 

1 (Xa — Yo.) (xp — yp) 

Ix—Ylop = Ix —y| bop eB B ’ 
(In|x—yl).0 = 72, 
|x—y/? 


(Xa —Ya) (xg —yp) 
|x—yl* 


1 
(In|x—y]) 0p = Ix — yl? ap 2 


? 


(Xa — Ya) (xp — yp) ; 
|x—y/? 


(lx—yl?In|x—yl) op = (2In|x—y| +1) bop +2 
consequently, from (2.7) and (2.12) it follows that 


(x1-y1)* 81 —y1) (%2 — ya) 


Al B 
puyel Pe Ivy? 
; _ _ _ _ 2 
cM yi) (x2 y2) Ali Ix y| pe y2) 
|x— yl? |x— yl? 


+O(|x—yl?In|x—yl), (2.13) 
where C is another combination of A and U. This indicates that for x close to y, 
D(x,y) = O(In|x—y|). (2.14) 


A similar computation, but with differentiation in terms of the point y, leads to 
the expansions 
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yi-*1 (vi —x1)(y2 — x2)? 
Pui (x,y) = vi(y) (an ae +b p—x4 
ad Oi ui)? — Oi a1) Ope) i 
y—aP | 
Y2— *2 (1 — x1)? (y2 — x2) 
+nby(apaet yal 
fy OTH) 2 2) — On 92)” — 
ly—-\4 ee 
2 
y2—%2 (v1 — x1)“ (2 = 2) 
Palsy) =") (an? bi2 pox 
We: x2)(yi — x1)? — (2 x) 
y— x14 
Yi=*1 (v1 = 1) 2-2)? 
+26) (diag tere 
(v1 = 1)? = (v1 = 1) (2 = 2)? 
fiz yal om 
2 — x2 (vi —x1)?(2 — x2) 
Palsy) =") (aor 235 bp) 
46, x2)(yi — x1)? = (2 Ba)? ses 
yaa? 
Yi=*1 (1 = 1) (2 — 2)? 
+b) (der 4 ba 
(1-1)? — G1 — x1) (2 — 92)? 
fr yal sai 
Yi-*1 (yi — x1) (2 — 2)? 
Pralssy) = ")) (aan —9y bpp UO, 
i) xi) — (v1 — 1) (v2 ty ats 
ly — x)? 
y2— X2 (y1 —x1)* (2 — 2) 
v2) (drape tea 
(1 — 1)’ (2 — 22) — 2-2)? 
+ fro aa ae 


where a, b,c, and d with various subs 
expansions show that for x close to y, 


P(x,y) = 


(2.15) 


cripts are combinations of A and w. These 


O(|x—y|7'). (2.16) 
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On the other hand, for y on 0S and |x| — ©, we have [6] 
Ko(h|x—y|) = O(|x|7/2e!) Va = const, 
which, by (2.12), (2.7), and (2.10), yields the asymptotic relations 
t(x,y) = O(|x[""7eM) 
D(x,y) = O(|x|7>/2e7 Fl) (2.17) 


The next assertion is the analog of Green’s representation formula. 


2.2 Theorem (Somigliana formulas). (i) Jf u € C?(S*) AC! (S*) and (Zu) (x) =0, 
x€ St, then 

u(x), xeESt, 

su(x), xe dS, (2.18) 
0, xeS. 


(ii) Ifu € C?(S~)NC!(S-) No and (Zu) (x) =0, x € S~, then 


[Penmno) — P(x, y)u(y)|ds(y) 


os 


0, xeSt, 
= [eenTo) —P(x,y)u(y)|ds(y) = 4 g(x), x€ dS, (2.19) 
as u(x), «ES. 


Proof. (i) Consider the region $+ \ Oye, Where Oy¢ is a disc centered at x and of 
radius €, and such that 6y,¢ C S* (see Fig. 2.2). 


x2 


x] 


Fig. 2.2 The domain St \ Gy¢. 
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Using Theorem 1.8 in S* \ 6,,¢ with the position of u and v reversed and v = v(y) 
replaced by the columns D() (y,x) of D(y,x), we find that 


[D'%)" (y, x) (Zu)(y) — ul (y)Z(%)D™ (y,x)] daly) 


S*\ One 


St\ Ove 


and the above equation becomes 


+ f DO" G.2)(Ta)Q) a7 9)T(G)D(.)]ds(9). 2.20) 


Taking (2.9) and (2.8) into account, we can write the integrands in (2.20) as 
Di)" (y,x) (Tu)(y) — u(y) (AD (9,2) 
x)(Tu)p (y) — ug (y)(T(A)D™ (y,x))p 
(Tu)g(y) — Tpy(Ay)Dya(y,x) up (y) 
(Tu)g(y) 


os 
+ f (D(ey)(Tw)(y)—PEy)uly)] 459) =0. 2.21) 
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Since u € C?(S+) C!(S*), it follows that all the partial derivatives of u are 
bounded on 0S, which means that Tu is bounded; therefore, by (2.14), 


/ D(x,y)(Tu)(y)ds(y) = O(eIne) 30 ase 0. (2.22) 


00x, 


Next, 


J Peex)ui)as() 
O0x.¢ 
= f PCy)uo)-ue dso) +4) f Peey)ds). 2.23) 


O0x.¢ I0x.¢ 


Going over to polar coordinates with the pole at x, we represent the point y € d0y.¢ 
as y = (€, @) and, taking into account the fact that the outward normal on 00,.¢ is 
directed into 0.¢, find that (2.16) has the more explicit form 


1 1 
P = 14+-¥¢(0 
(y) = 7 I+ 5 98), 
where 


2n 
[9(e)a0 ~0. 
0 


Then, for y € 00x.¢, 


P(x,y)=O(e'), u(x)—u(y) = Ole), ds(y) = ed, 


/ P(x,y)[u(y) —u(x)]ds(y) = Ole) O(e)e 30 as e 0. 


00x O0x,€ 
20 
= u(x) (1+ [-9(0)40) = u(x) 
0 
Substituting this in (2.23), we conclude that 
/ P(x,y)u(y) ds(y) > u(x) ase > 0. (2.24) 


O0x.€ 
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The first part of (2.18) is now obtained by replacing (2.22) and (2.24) in (2.21) and 
letting € > 0. 

If x € OS, only one half of the disc 6,.¢ is needed to isolate the point x. If x € S~, 
D(x,y) is analytic at any point y € St, and the argument involving the disc 0;.¢ is 
not needed. This generates the second and third parts of (2.18). 


(ii) Let Kr be a disc with the center at x and radius R sufficiently large so that 
S* C Kp, and let S* = Kp \ S* (see Fig. 2.3). Using the first part of formula (2.18) 
in S* = S* VOSUOKg and noting that in this case the outward normal on 0S points 
into St, we have 


— [[D(,y)(TH)(y) — Pla y)u(y)] ds(y) 
os 
+ f Dey)(Tw))— Pley)uy)]ds(y) = ul). 
OKR 


In view of the asymptotic relations (2.17), the second integrand vanishes as R — 9, 
and we obtain the third part of (2.19). The argument for x € 0S and x € S* is similar 
to that in (i). 


x] 


Fig. 2.3. The domain S*. 
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2.3 Definition. The functions 


Voy) = [ Dixy)@0) avy), 
os 


(W(x) = [Peer wO)ds), 
os 
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where @ and y are 2-component vector functions defined on OS, are called the 
single-layer potential and double-layer potential of density @ and yw, respectively. 


2.4 Theorem. /f 9, wy € C(0S), thenVo,Wy Ee &. 


Proof. The asymptotic behavior (1.12) follows immediately from (2.17). 


2.5 Theorem. (i) Jf @, yw € C(AS), then Vp and Wy are analytic in S* US~ and 
Z(Vo)=Z(Ww)=0 inStUS-. 
(ii) If 9 € C° (AS), then Vp € C-*(R?). 
(iii) If y € C°-* (0S), then W*y € C°(S*), where 
Wry = (Ww)|se. 


(iv) If @ € C°* (0S), then the direct values Vo~ and Woy of Vp and Wy on 0S 
exist, with the integral in Woy understood as Cauchy principal value. 


(v) The operators ¥* and ¥~ defined by 


V9 =(Vo)|5+, 
“p (Vp) |s+ (2.25) 
V9 = (V@)|5- 
map C°-* (QS) to C= (S*), a € (0,1), respectively, and 
T(V*o) = (We +i De, 
(Y¥*9) = (We +351) (2.26) 


T(V~ 9) = (Wo -31)@, 
where W,, the adjoint of the operator Wo, is defined by 
(W5 o)(x) = [(TA)D(.»)) OO) ds(y),_-xE aS, 
os 
with the integral understood as Cauchy principal value. 
(vi) The operators W* and W~ defined by 
yy LW Vlss ins, 
(Wo—5l)w onds, 
( 


7 Ww)ls- in S-, 
W-w= 
. eee on 0S 


map C°-* (QS) to C= (S*) and C!}“ (AS) to C\™(S*), a € (0,1), and 
TWrw)=TW-w) onc*(dS). 


(vii) The operator Wo maps C°*(0S) to C!™ (AS), a € (0,1). 


20 2 The Layer Potentials 


Proof. Statement (i) is verified by means of Theorem 2.1. 

Comparing the asymptotic expansions (2.13) and (2.15) with the corresponding 
ones in the absence of an elastic foundation [6], we see that they contain the same 
terms but with different coefficients. Since properties (ii)—(vil) depend exclusively 
on the order of the singularities in D(x,y) and P(x, y) for x close to y, their proof is 
a mirror image of that of their equivalents in [6]. 


2.6 Remarks. (i) In Theorem 2.4, the derivatives on 0S of the functions defined on 
S+ and S~ are one-sided. 


(ii) Clearly, 
Z(V~ =0, xeS', 
(Y*y(x)=0, x oo 
ZV p)(x) =0, xES. 
(iii) If 
Vi0=V o=), 
then 
o=0. (2.28) 


This is easily seen from the fact that 
T(V")p=T(V)p=0, 


and (2.28) follows from (2.26). 


(iv) The Somigliana relations can be written in a more compact form in terms of 
the operators /* and W*. Thus, if u is a solution of Zu = 0 in ST, then 


V*(Tu)—-W* (ulas) =u, V~ (Tu) —-W~ (ulas) = 0, (2.29) 
and if wis a solution of Zu =0in S-, then 
—V~(Tu)+W (ulas) =u, —V¥*(Tu) +H (ulgs) =0. (2.30) 
(v) From Theorem 2.5(vi) it follows that we can define a boundary operator 
No: Ch%(0S) + C°*(dS) 


by setting 
Now=TWy)=TW y), ywecr"(ds). (2.31) 


(vi) The operators defined in Theorem 2.5(v),(vi) satisfy 


(Y @)las=% P =Vo9, 
(W*Wlas = %o W = (Wo# 5!) V. 


(vii) We note that Vo, Wo, and W; are defined on C°:“(S), whereas No is defined 
on Cl (QS), a € (0,1). 


(2.32) 
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2.3 Properties of the Boundary Operators 


The operators defined in the previous section have a number of important properties 
that will be used extensively in the development of solutions of our boundary value 
problems. 


2.7 Theorem. Vo, Wo, Wo , and No satisfy the composition formulas 


WoVo =VoWo, NoVo=Wo?-4G1 onC®*(dS), (2.33) 
NoWo =WeNo, VoNo=Wo—41 onc!*(as). (2.34) 

Proof. Let 
u=Vto-Wy, (2.35) 


where @ € C° (0S) and y € C!(QS) are arbitrary. By Theorem 2.5(i), (Zu) (x) =0 
for x € ST. Using (2.32) and restricting u to 0S, we get 


ulas =(V* @)las— (WT Was; 


or, setting ulas = a, 
a =Vop — (Wo —51)y. (2.36) 


Applying the boundary stress operator T in (2.35) and using (2.26) and (2.31), we 
arrive at 


Tu=T(V*o)-T(Wtw); 


that is, with Tu = B, 
B = (Wo —51)0—Noy. (2.37) 


Next, we write the Somigliana representation formula (2.29) restricted to the 
boundary OS: 


(V*(Tu))|as— (W* (ulas) las = 4las: 
which is the same as 
Vo(Tu) — (Wo — 3 1)ulas =ulas, 


or 
VoB — (Wo- 5 1a=a. (2.38) 


Finally, applying T to (2.29) and using (2.26) and (2.31), we have 
T(¥* (Tu))— TW" (ulas)) =Tu, 


or 
(Wi +31)B—Noa =B. (2.39) 


System (2.36), (2.37) can be written as 
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(5! tn) (@)-C) am 


Also, system (2.38), (2.39) takes the form 


(5! t)@)-@) om 


Replacing (8, a)! from (2.40) in (2.41), we now deduce that 


(ines —No ) (2) = ("oH —No )(2) 
Yo —(Wo-51)/) \v Yo —(Wo-32)/) \w/)’ 


and, since the matrix on the left-hand side expands to 


( (Wox+h1)?—NoVo — —(We +42) No+N0(Wo _ 
Vo(Wg +41) —(Wo-41)Yo — —VoNo-+(Wo— 41)” | 


the usual identification of elements and the arbitrariness of @ and y in their corre- 
sponding spaces lead to (2.33) and (2.34). 


2.8 Theorem. Wo + 51 and Wy + 51 are operators of index zero. 


The proof of this theorem is very similar to that of Theorem 2.5(v) and of the 
corresponding assertion in [6], where the terms with the strongest singularity in the 
kernels of Wo and W, are the same as here, up to a constant factor. Since the index [6] 
is homogeneous of order 0 in the numerical coefficient of 0,1n |x — y|, we conclude 
that the index of our operators Wo and W) is the same as that of their counterparts 
in [6], namely 0. 

A useful tool from functional analysis that will enable us to solve our boundary 
value problems is the Fredholm alternative, which can be stated using Definitions 
1.72 and 1.73 from [7]. 


2.9 Definition. A dual system (X ,Y) is a pair of normed spaces X and Y with a 
non-degenerate bilinear form (-,-):X x Y >C. 


2.10 Definition. Let (X,Y) be a dual system with bilinear form (-,-), let K :X > X 
be an operator with a (unique) adjoint K* : Y > Y, let@ € C, m £0, and consider 
the equations 


(K*-@l)p=g, geEX, (K*) 


together with their homogeneous versions (Ko) and (Kj). We say that the Fredholm 
alternative holds for these equations in (X,Y) if either 


(i) (Ko) has only the zero solution, in which case so does (Kj), and (K) and (K*) 
have unique solutions for any f € X and any g € Y; or 
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(ii) (Ko) and (K) have the same finite number of linearly independent solutions 
{P1,21°-° Pn} and {Wi Yo,---, Wn}, 
respectively, and (K) and (K*) are solvable if and only if 
(f,Wi)=0 and (g,9,)=0, i=1,2,---,n. 

2.11 Corollary. The Fredholm alternative holds for the pairs of equations 

(Wo-3o=f, (Wo-Z)w=s, 

(Wot+sNe=f, (Wtz)w=s 
in the dual system (C°-“(0S),C°*(AS)), a € (0,1), with the bilinear form 

(PW) = J otwas. 
os 


2.12 Theorem. The null spaces of Wo — 51, Wi — 51, Wot 41, Wi +41, No, and 
Vo consist of the zero vector alone. 


Proof. If @ € C°-*(AS) is in the null space of We" — ae then, by (2.26), 
(Wo — 31) 9 =T(Y 9) =0 
and, by (2.27), Z(V~ @) = 0; therefore, 
V-o=0 


as the unique solution of the homogeneous problem (N7). By (2.25) and (2.32), 
Y~ o = 0 implies that 

(Y* @)las=(% 9)las =9 
and, by (2.27), Z(V* @) =0. Consequently, Y * @ = 0 as the unique solution of the 
homogeneous problem (D*). Since 


V+o=0, V-o=0, 


from Remark 2.6(iii) it follows that @ = 0. 
If p € C(AS) is in the null space of W; + 57, then, by (2.26), 


(Wo +31)9 =0 


is equivalent to T(¥*@) = 0 and, by (2.27), Z(V¥*@) = 0; therefore, Yo = 0 
as the unique solution of the homogeneous problem (N*). By (2.25) and (2.32), 
¥* p =0 implies that 

(Y ®)las=(V%* ®)las =0 
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and, by (2.27), Z(V~ @) = 0. Therefore, since 
(VY ®)las=0, V QE, 
¥Y~ o = 0 as the unique solution of the homogeneous problem (D7). Given that 
V*0=V-o=0, 


from Remark 2.6(iii) it follows that @ = 0. 
Corollary 2.11 and the Fredholm alternative imply that the null space of Wo + 51 


and that of Wo — 51 also consist only of the zero vector. 
If Nop = 0, then, by (2.34), 


0 = Vo(No@) = (VoNo) 
= (Wi — 41) 9 = (Wo— 31) (Wot 31)9, 
so, given the structure of the null spaces of Wo — 51 and Wo + 51 , we conclude that 
o=0. 
By (2.32), Vop = 0 implies that 


(Y* @)las =0. 


By (2.27), Z(V + @) =0; hence, by Theorem 1.9, “+ @ = 0 as the unique solution 
of the homogeneous problem (D*). Also, by (2.32), Vo@ = 0 implies that 


(Y )las =9. 
By (2.27), Z(V~ @) = 0; therefore, by Theorem 1.9, since 
(VY @)las=0, V DEW, 


Y~ o =0 as the unique solution of the homogeneous problem (D~ ). The equalities 
V+ op =V~@ =0 and Remark 2.6(iii) now imply that ~ = 0. 


Chapter 3 
Existence of Solutions 


In this chapter, we describe and apply two boundary integral equation methods to 
solve the fundamental boundary value problems formulated in Sect. 1.2. 


3.1 The Classical Indirect Method 


For the interior and exterior Dirichlet, Neumann, and Robin problems (D~), (N*), 
and (R*), the solutions are sought in the form 


u=Wr@ for (D*), 
u=W for (D-), 
u=Vw for (N*), 
u=V-w_ for (N-), 
u=VTo for (R*), 


u=¥~@ _ for (R-). 


From the properties of W* and ¥Y* listed in Theorem 2.5 and the boundary con- 
ditions of each of these problems, the unknown density functions @ and y must 
satisfy, respectively, the boundary integral equations 


(Wo-3No=F, (Qe) 

(Wo + 4 No=&, (Qe) 

(Wo +31) w= 2, (4°) 

(Wo -3)y=F, (4) 

(Wo +31) 9+0(Vo@) = (Zé) 

(Wo — 31) 9—6(Vo@) = (Zo) 
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3.1.1 The Dirichlet Problems 


3.1 Theorem. (i) (Zé) has a unique solution p € C!*(S) for any P € C'*(dS), 
a € (0,1). Then (DT) has the unique solution 


u=W9@. (3.1) 


(ii) (Fo) has a unique solution g € C!*(AS) for any Z € C! (dS), a € (0,1). 
Then (D7) has the unique solution 


u=W~9. (3.2) 


Proof. (i) By Theorem 2.12, the null space of Wo — 51 contains only the zero vector, 
so, by the Fredholm alternative, (Yt) has a unique solution p € C! (dS). 
Since u defined by (3.1) satisfies Zu = 0 and 


ulas = Wir 9 = (WoL =P, 


we conclude that this function is the unique solution of (D*). 

(ii) The proof of this assertion is a mirror image of that of (i). Additionally, here 
we must mention that, by Theorem 2.4, u defined by (3.2) also satisfies the far-field 
condition u € &. 


3.1.2 The Neumann Problems 


3.2 Theorem. (i) (4%) has a unique solution y € C! (0S) for any 2 € C°“ (dS), 
a € (0,1). Then (NT) has the unique solution 
u=Vty. (3.3) 


(ii) (WC) has a unique solution y € C° (0S) for any S € C° (dS), & € (0,1). 
Then (N— ) has the unique solution 


u=V-w. (3.4) 


Proof. (i) By Theorem 2.12, the null space of Wy + 51 contains only the zero vector, 
so the Fredholm alternative implies that (.4.*) has a unique solution y € C'(dS). 
Since u defined by (3.3) satisfies Zu = 0 and 


Tu=T(¥ ty) = (Wi +43)W=2, 


we conclude that this function is the unique solution of (NT). 


(ii) The proof of this part is similar to that of (1), with the added mention that, 
according to Theorem 2.4, u defined by (3.4) also satisfies u € &. 
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3.1.3 The Robin Problems 


3.3 Theorem. (i) (Z<) has a unique solution p € C°*(0S) for any H € C° (dS) 
and any 6 € C°*(0S), a € (0,1), Then (Rt) has the unique solution 


u=V*9. (3.5) 


(ii) (@e) has a unique solution p € C°*(AS) for any Z € C°“(AS) and any 
oEeC® (95), 0 € (0,1). Then (R7) has the unique solution 


u=V' 9. (3.6) 
Proof. (i) Let p be a solution of the homogeneous equation (Be ); that is, 
(Wo +32) H+ 0(Vou) =0 
which, according to Theorem 2.5(v), means that 
TV p)+o(%h) = 


Hence, u = Y* p is a solution of the homogeneous problem (R*), so, by Theorem 
1.9, ¥*p =0. Then, by Theorem 2.5(ii), 


% H=H' n=O. 


Since Z(V~ 1) = 0 and, by Theorem 2.4, Yu € &, we conclude that V~ p is 
a solution of (R~), so, by Theorem 1.9, Y~ u = 0. Remark 2.6(iii) now implies 
that u = 0. Consequently, by the Fredholm alternative, (Ze ) has a unique solution 
@ €C°(9S). 

Given that u defined by (3.5) satisfies Zu = 0 and 


Tut oulas =T(V +9) +0(%'o) = (We +5 e+0(Ho) =~%, 


this function is the unique solution of (R*). 


(ii) The proof of this assertion is similar to that of (i), with the added observation 
that u defined by (3.6) also satisfies the far-field condition u € .&. 


3.2 The Direct Method 


In this method, we use the Somigliana representation formulas to set up appropriate 
boundary integral equations for (D~), (N~), and (R*). 
For the interior problems, the restriction of the first equation (2.39) to 0S yields 


Vo(Tu) — (Wo + 51) (ulas) = 0. 
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In (D*), ulas = F is known and Tu = @ is unknown; in (N*), Tu = J is known 
and u|j; = y is unknown. Hence, the corresponding boundary integral equations 
are 


Yop = (Wot 5D F, (Fp) 
(Wo +431) y=VWo2. (457) 


For the exterior problems, the first equation (2.40) restricted to 0S is written as 
—Vo(Tu) + (Wo — 1) (ulas) =0. 


In (D~), ulas = Z& is known and Tu = @ is unknown; in (N_), Tu = -Y is known 
and u|j9; = y is unknown. Hence, the corresponding boundary integral equations 
are 


Yop = (Wo— 3&2, (Dp) 
(Wo 31) y=WF, Crs 


Both (Ys ) and (A, ) are equations of the first kind, whereas (.4%") and (.% ) 
are equations of the second kind. 
For (R*), we start with the boundary condition written in the form 


Tu = —oO(ulg5)+-%, 
which, on substitution in the first equation (2.29), leads to 
u=V*(—o(ulas) +H) — W* (ulas). 
Restricting this equation to 0S, setting 
y =ulas, 
and using (2.32), we find that 
y=—Vo(ow)+ (Wo- sl) yw-V(%), 


or 
Vo(ow) + (Wo+ sl) y=Vo%. (2p) 


For (R7), writing the boundary condition in the form 
Tu =O(ulas) +2 
and using a similar argument, we arrive that 


Vo(ow) — (Wo- 41) y=-WZ. (Zp) 
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3.2.1 The Dirichlet Problems 


3.4 Theorem. (i) (G3) has a unique solution g € C* (AS) for any P € C!(dS), 
a € (0,1). Then (DT) has the unique solution 


u=Vt9-—W*®. (3.7) 


(ii) (BH) has a unique solution g € C°“ (AS) for any R € C!™ (AS), « € (0,1). 
Then (D7) has the unique solution 


u=—-V Qo+W &. (3.8) 
Proof. Applying No to (Ys), we get 
No(Vo~) =No((Wo +51) P). 
By Theorem 2.7, 
No(Vop) = (Wo? — £1) @ 


and 
No| (Wo + 51) P|] =W*(NoP) + 3NoFP = (W* +51) (NP); 
hence, 
(Wo? -i Do = (W* +41) (NP), 
or 
(Wi +50) [(Ws - 3D e-NF] =0. 


By Theorem 2.12, the null space of Wy + sl consists only of the zero vector; there- 
fore, 
(Wo — 31) 9-NoP =0. 


Since, by the same theorem, the null space of Wy — 51 also consists only of the zero 
vector, from the Fredholm alternative it follows that the above equation has a unique 
solution @, which is also the unique solution of (DB ). 

Given that u defined by (3.7) satisfies 


Zu=0 inSt 
and, in view of (Fz), 
ulas = (VT @)|as— (W* Plas 
=Vop — (Wo — 51) P 


we conclude that (3.7) is the unique solution of (D*). 


(ii) The proof of this assertion is similar to that of (i), with the additional mention 
that, by Theorem 2.4, u defined by (3.8) also satisfies u € 2. 
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3.2.2 The Neumann Problems 


: 5 Theorem. (i) (-4%" ) has a unique solution y € C!\“(0S) for any 2 € C°* (0S), 
€ (0,1). Then (Nt) has the unique solution 


u=Vt*Q-Wy. (3.9) 


(ii) (4%) has a unique solution y € C! (0S) for any S € C° (dS), @ € (0,1). 
Then (N— ) has the unique solution 


u=-V~Q4W-w. (3.10) 


Proof. (i) By Theorem 2.12, the null space of Wo + 51 consists only of the zero 
vector, so, by the Fredholm alternative, (./%;") has a unique solution y € C!'“(0S). 


Clearly, u defined by (3.9) satisfies 


Zu=0 inS*. 


To show that it also satisfies the boundary condition Tu = 2, we apply T to both 
sides in (3.9): 
Tu = (Wo +31) 2—Noy. (3.11) 


Let 
H = (Wo +51 2-Noy— 2= (Wy — 51) 2-Now. 


Operating on both sides above with Vo and using (2.34), we see that 
Vo H = Vol (Wo — 51) 2-Noy] 
= (VoWj )2—-4V2— (Wi - 41) 
= (WoY)2-5V2-(We- FI 
= (Wo $1)2— (Wey 
= (Wo — 51) [Yo2— (Wo- 51) y] =0. 
Therefore, by Theorem 2.12, # = 0, which leads to 
(Wy +51) 2-Mw-2=0, 


or 
(Wo +5 2-Ny=2; 
hence, in view of (3.11), 


Tu= 2. 


(ii) The proof of this part is similar to that of part (i), with the additional obser- 
vation that, by Theorem 2.4, u € &. 
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3.2.3 The Robin Problems 


3.6 Theorem. (i) (#5) has a unique solution y € C!\“(0S) for any H € C°“ (dS) 
and any 6 € C°*(0S), a € (0,1). Then (Rt) has the unique solution 


u=—Vt(ow)-Wtwtv' #. (3.12) 


(ii) (5) has a unique solution w € C!\*(AS) for any LZ € C°“(dS) and any 
o €C**(9S), a € (0,1). Then (R~) has the unique solution 


u=—-V (oW)+W w-V X. (3.13) 
Proof. (i) The homogeneous version of (4%, ), namely 
Vo(ow) + (Wot+5l)w=0, 
can be rewritten as 


Vo(oy) + (Wo- 51) w= —v. (3.14) 


Let 
U-=V (ow)+W (y). 


Then ZY~ =OinS” and 
WU \as =Vo(oy) + (Wo+5l)y, 


which means that Y~ is a solution of the homogeneous problem (D~ ), so, by The- 
orem 1.9, 


uU =0. 
Next, applying T to Y~ and using (2.26) and (2.32), we obtain 
T(Y (ow) +T HY (y)) =, 


so 
(Wi — 41) (oy) +Now=0. 


or 
(Wy +51) (oy) +Noy = oy. (3.15) 


Let 
U*=V"(ow)+W* (yp). 


Then ZY* = 0 in S* and, by (2.26), (2.31), (3.15), and (3.14), 
TU* +0U* \g9=T(Vt (oy) +Wty| + o[¥* (ow) + ¥*(W)llas 
= (Wo +31) (oy) +My +o[Vo(oy) + (Wo- 51) y] 
=oy—oy=0. 
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Hence, Y%* is a solution to the homogeneous problem (R*), and from Theorem 1.9 
it follows that Y@+ = 0, which means that the homogeneous equation (Be ) has a 
unique solution. Then, by the Fredholm alternative, so does (2B ). 

It is obvious that u defined by (3.12) satisfies 


Zu=0 inS?t. 
Restricting (3.12) to AS and using (4; ), (2.26), (2.31), and (2.32), we see that 
Tu = —(Wo +31) (oy) —Now + (Wo + 3D) # 
= (We + 51) (4% —oy) —Noy, 


so 
Tu+ oulas = (Wo +51)(% —oy) —Noy 
+o6[Vo(—oy) — (Wo — 51) w+Y.#] 
= (Wo + 51) (4% —ow) —Now 
+0[Vo(.% —oy) — (Wo + sl)w+ y| 
= (Wo +431)(%# —oy)—My+ov. (3.16) 
Let 


PY = (Wo +31)(# —oy)—Noy+oy-# 
= (Wp — 31)(% —oy)—Noy. (3.17) 
By (#5), (2.33), and (2.34), 


Vo = Vo[(Wo — 31) (4% — oy) —Noy] 
= Vo(Wo — 31) (4 — ow) —Vo(Now) 
= (Wo — 31)Vo( H — ow) — (We - 5D 
= (Wo — 31) [Vo(-% — ow) — (Wo+ 51) yw] =0. 


According to (2.32) and Theorem 2.12, Vo'¥ = 0 implies that Y = 0; therefore, by 
(3.17) and (3.16), 


(W5 +51)(% —oy)-Myt+ov=%, 
or 
Tu+oul\j5 = 4%, 


which confirms that (3.12) is the unique solution of (R*). 


(ii) The proof of this part is similar to that of part (i), with the additional detail 
that, by Theorem 2.4, the function defined by (3.13) also satisfies u € &/. 
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3.7 Remark. An alternative way to reduce (R*) to a boundary integral equation is 
to write the boundary condition in the form 


ulas = o | (XH — Tu), 


which, replaced in the first equality (2.29) restricted to 0S and combined with (2.32) 
and the notation 
p =Tu, 


leads to 
Yop + (Wot 31)(o '@) = (Wot 5I)(o |X). (Bx) 


A similar equation can be set up for the problem (R7), and an existence theorem 
analogous to Theorem 3.6 can be proved. 


3.8 Remarks. (i) The existence and uniqueness results established in this chapter 
remain valid if the boundary 0S is closed, continuous, and piecewise smooth—for 
simplicity of the analytic argument, in our case we assume piecewise continuous 
tangent and curvature—since this allows us to apply the divergence theorem. How- 
ever, in such domains we expect the solution to lose some degree of smoothness. In 
a rectangle, for example, this happens as we approach a corner, even if the Dirichlet 
data are continuous there. Chapter 5 contains several numerical approximations for 
domains of this type. 


(ii) The direct and indirect methods can also be adapted to the case when the 
boundary conditions are of mixed type—for example, when Dirichlet conditions are 
prescribed on one part of the boundary and Neumann conditions on the remaining 
part. 


Chapter 4 
Software Development 


4.1 Programming Environment 


All computational, symbolic, numerical, and graphic software used in the rest of the 
book has been generated with Mathematica®, a commercial product of Wolfram 
Research. The purpose of this chapter is to both document and illustrate how the 
software has been developed. Many issues specific to the boundary integral meth- 
ods that form the object of our study require special attention. Below, we summarize 
these issues and indicate how certain difficulties associated with them can be over- 
come in the Mathematica® environment. We do not intend to give specific detailed 
code, but rather a general outline of how the code has been assembled. 

To simplify the terminology, in what follows the words ‘code’ and ‘coding’ refer 
to programming with Mathematica®. Also, the indentation of the displayed formu- 
las has been changed to conform to the TEX output of the software. 


4.2 Notation 


Mathematica® is at its core a functional language. All functions and their arguments 
are given a name coded in the format 


namelarg1,arg2,...,argn]. 


This allows great flexibility in generating functions that accomplish the intended 
tasks. However, such notation is both cumbersome and, often, difficult to follow, 
which has prompted us to replace the most commonly used functional calls in our 
software by suggestive symbolic notation. 


4.1 Example. It is very convenient to code Mathematica® definitions for the usual 
infix operators as x, 0, and ©, respectively, to denote the scalar, functional, and 
operator compositional ‘inner products’. Thus, the composition of a matrix of oper- 
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ators with a matrix of functions by means of the operator © is written as 


isl#l& firle]&) | failel\ _ f fialsilell + falsolel] 

Ar ifA& fo2l#]& g2I¢] failgilél] + Algol] 
Here, the functions are entered in Mathematica® notation as pure functions by 
means of the customary notation 


fi jl#1&. 


This type of symbolic notational definitions can be combined with others to build 
complicated expressions that take on a more traditional mathematical appearance. 


4.2 Example. If for the line integral we introduce 


Notation] f f_dij_ << __ linelntegrall[f_,t]], 
r 


then in our code we can use the familiar form 


¢ fltlar 
Tr 


instead of its rather obscure equivalent 
lineIntegralI’ [f(t], 1]. 


4.3 Example. The representation formula (2.18) for x € S* quoted later in this chap- 
ter is coded as 


ulx] = f Dbx velo (Tu)(yll] x dsat( ar 
Tr 
ee $ Plx,y[t]] ouly[e]] x dsdt{¢] any. 
Tr 


In what follows, expressions related to code in Mathematica® are displayed as 
much as possible in this suggestive notation. Also, whenever appropriate, specific 
functional calls will reference equations and formulas developed earlier in the book, 
to aid the readers in their interpretation of certain important details. 


4.4 Remark. The text in Chaps. 4 and 5 contains a mixture of mathematical sym- 
bols and Mathematica® objects. To distinguish between the two, the former are 
written in the standard mathematical font used in the first three chapters, and the 
latter in a combination of fonts according to the conventions embedded in the soft- 
ware. Sometimes, both categories of symbols need to stand side by side in the same 
formula. 
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4.3 Coding the Mathematical Model 


The system of equations describing the equilibrium state of a thin plate on an elastic 

foundation was described in Chap. 1. These equations, defined in regions S* and 

S~ with boundary 0S (see Fig. 1.1), form the linear and homogeneous elliptic sys- 

tem (1.8). In Mathematica®, using notational definitions, we enter the differential 

operator (1.6) as 

_ M(x, xy # + Oxy xp #) + (A + Mh) Oey xy KAR (A + BH) Ox xp #& 
© \ (A+ Hr antt& H(A, 01# + Ora xa) + (A+ H) Orn writ — AB] ” 

(4.1) 


where the subscript x on the left-hand side indicates that the operator Z acts with 
respect to the point x. This subscript becomes relevant when Z is applied to two- 


point functions. 
Similarly, the coding of the boundary stress operator (1.10) is 


_ (A +2m) V1 [x1,x2]0x,#+ UV2[1,x2],#& = V2[x1,x2]0r, #+AVI [x1 , x2], #& 
“VN Av2 [a1 x2], H+ MV fr x2] A #& UV A [x12], #+ (A +2) v2 [ar 2] #& J 
(4.2) 


The homogeneous system (1.8), now rewritten as 
Z,.04=0, xES™US, 


together with Dirichlet, Neumann, or Robin boundary conditions prescribed on 0S 
constitute a boundary value problem (see Sect. 1.2). In the rest of the book, we focus 
our attention exclusively on problems in S* US and, for brevity, denote St by S. 


4.4 Coding of the Matrix of Fundamental Solutions 


The symbolic capabilities of Mathematica® allow us to compute the explicit form 
of the matrix of fundamental solutions constructed in Sect. 2.1. As shown below, 
these solutions are nontrivial, and the software both removes the potential for deriva- 
tion errors and makes it possible for us to verify the result. In our case, the 2 x 2 
matrix of fundamental solutions D satisfies (see (2.1)) 


= Di [x,y] Dizlx,y]\ 10 
Z,° D|x,y] = Zz 0 Be aad = —6[x,y] iG . . 


For convenience, we use our own coding definitions to recall the procedure for 
constructing D. According to (2.2), D is of the form 


D{x,y] = adjoint{Z,] (’ >I me) (4.3) 
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which, substituted in the equilibrium system, yields 


Zy © adjoint|Z,] (’ fa ‘ - = Della? (’ . i) 


ae er (; ae 


The matrix of fundamental solutions is therefore constructed from the solution of 
the equation 


Det[Z;] oz[x,y] = —6 [x,y]. (4.4) 
Mathematica® evaluates the symbolic determinant operator Det [Z,] as 
(K7#1 —K(A +31) Ox] xp #1 — (A+ 3b) Ox yo#l + AMO x1 xd xt #l 


2 2 
+2u Ont xi xixitl + 2A MOx] x1 x2x2#1 +4u Ox] x1. x2.x2#l 
2 
Tv AMOx? x2.x2.x2#1 +2u x2 x2.x2.x2#1) &, 


which simplifies to 
(PI—KA+3p)A+U(A+2p)AA)&; (4.5) 


here, J denotes the identity operator/matrix since the symbol /, used in the earlier 
chapters, has a different assignment in Mathematica®. 
This operator can be formally factored, so the equation for t becomes 


k k 
A+2u)(A T){A 3) ot[x,y] = —6[x,y]. 
wa +2u)(4-£3) (4-7 523) onfs9] = Stn 
Converting the above equation to polar form, solving it symbolically, and keep- 
ing, as indicated in Sect. 2.1, only the terms with singularities at r = 0 leads to ¢ 


being expressed as the combination of modified Bessel functions (see (2.5)) 


1 i : 
t{r] = Teak LH) (Besseik o nf] — BesselK 0 va | ) . 


At this time, our code can be used to verify the accuracy of this expression by 
checking symbolically that (4.4) is satisfied. 
All the necessary elements are now available to compute the matrix of funda- 
mental solutions (4.3), where (see (2.6)) 
adjoint [Z,] 
—k#1+ Hox] xy #l +(A +2) 0x9 yp#1& —(A + Hx] x2#1& 
~ Gey KH + (A+2u) Oxy x #1 ivan 
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Given the excessive size of the full expression, we show explicitly only four of the 
60 terms that make up the component Dj of D: 


D,1[x1,x2,yl,y2] =--- 
BesselK [o. Vly? + (2 = 2/4] 
2n(A +H) 
x2y2uBesselK |0, /(x1 yl)? + (x2 ae 
m((x1 —y1)* + (x2—y2)?) (A +p)(A +2p) 
2x2y2u y/ cA Besselk |1, (xl — yl)? + (X2—y2)?\/74o| 
kn ((xl — yl)? + (x2—y2)?)°/? (a +n) 
x2y2ABesselK [2, /(x1— yl)" + (X2—y2)/4] 
2m ((x1 — yl)? + (x2—y2)?)w(A +H) 


It can be partially verified that this expression is correctly computed by using the 
symbolic capabilities of Mathematica® to check that 


Zxo D[x,y] = ¢ a) , xy. (4.6) 


As expected, D;,; and D2 have logarithmic singularities. Mathematica® con- 
firms this: when instructed to expand each of the four components, converted to 
polar form, in a Taylor series (see (2.13)), the result indicates that 


_ (OlLog{r]) oft] 
nul Gi oh eetal) | 4-7) 


4.5 Remark. Our boundary integrals are of the form 


$ Kix, yifblasbl, (4.8) 
Os 


where the elements of the matrix kernel K may have finite jump discontinuities or a 
weak or strong singularity at x = y. Theoretically, since y € 0S, 


(i) when x € S, (4.8) is a proper definite integral that can be computed without 
difficulty, however close x is to 0S; 

(ii) when x € OS, (4.8) is computed as an improper or a Cauchy principal value 
integral. 


The situation is less clear-cut if (4.8) is evaluated numerically. Mathematica® uses 
floating point arithmetic, so the coordinates of a point are always computed with a 
certain degree of inaccuracy. This means that a point x lying in S very close to the 
boundary may be placed numerically on 0S, or, the other way around, a point y on 
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OS may be located numerically inside S in the immediate vicinity of the boundary. 
Consequently, numerical evaluation of (4.8) may or may not conform to cases (i) 
and (ii) above. It is therefore necessary for us to examine the behavior of K[x, y] in 
all possible permutations: both x and y in S, both x and y on 0S, x in S and y on OS, 
and x on 0S and y in S. 


4.6 Example. Figure 4.1 illustrates graphically how each component of the ma- 
trix of fundamental solutions D[x, y] behaves if both x and y are in S. Taking, say, 
(x1, x2) = (0.5, 0.5) as a test point inside the domain, we notice that D, 1 [x,y] and 
D2 [x,y] have a logarithmic singularity and that Dy 2[x,y] and D2 [x,y] have a finite 
jump discontinuity at (yl, y2) = (0.5, 0.5), whose size depends on the direction of 
the path taken to approach the singular point. These graphs confirm our expectation 
that the numerical computation of boundary integrals with kernel D[x,y] will not 
encounter major difficulties when the singular point is in the interior of the domain. 
However, this is not so when the singularity lies on the boundary. 


Fig. 4.1 The components of D[x,y], x,y € S. 
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4.7 Example. Consider again the point (x1, x2) = (0.5, 0.5), but this time for a 
problem where the boundary is the line yl = 0.5. Figure 4.2 shows the behavior of 
the four components of D[x,y] when the point (yl, y2) is close to (0.5, 0.5). Their 
graphs are displayed for the domain to the left of the boundary and truncated to the 
right of it. The logarithmic singularity of D, 1 [x,y] and D2 2|x, y], evidenced by their 
graphs, becomes an important issue when we use numerical techniques to evalu- 
ate boundary integrals of the form (4.8) with K[x,y] = D|x,y]. Special care must be 
taken in such cases to preserve the accuracy of the computed result. We also see 
that the graphs of D, [x,y] and D21|x,y] have a finite jump discontinuity at x = y. 
Analytic arguments indicate that this does not happen if the boundary curve has a 
continuous tangent. For the specific case shown above, the corresponding graphs 
illustrate the lack of a jump discontinuity along the boundary curve at the point 
(yl,y2) = (0.5,0.5). However, unlike the purely analytic treatment, numerical inte- 
gration is not exact, and if the calculated point (yl, y2) is slightly off the boundary at 
y2 = 0.5, a jump discontinuity occurs, which can produce significant computational 
errors. 


ligt 
{] sl 
ONY) 


Fig. 4.2. The singularity in the components of D[x,y], x, y € OS. 
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Next, we construct symbolically the 2 x 2 matrix P, which is used to define the 
double-layer potential (see (2.9)): 


Pix, y]:=(ZoDpy,x])”. (4.9) 


We do this by direct substitution of the relevant quantities in the above formula, 
taking care of the necessary switch between x and y. The expression of P is very 
large, with a LeafCount of 29,183, which exceeds that for D by a factor of 10. 
In ordinary terms, this means that writing in full the four components of P would 
produce over 40 pages of printout. This underscores the importance of the symbolic 
computational abilities of Mathematica®. 

To better understand the internal structure of P, here we show only a small portion 
of P,.1; specifically, we list a sample of four of the 240 terms that make up this 
component: 


P,1[x1,x2,yl,y2] =--- 


3x1°y2uBesselK [o. J (xl = yl)? + (x2 —y2 2 /#] v2[yl,y2] 
4m ((x1 — yl)? + (x2—y2)2)” (A+) 
3xlyly2u? \/ i BesselK |1, /(x1 yl)? + (x2—y2 2a v2[yl,y2] 
ke (x1 — yl)? + (x2— y2)2)° (A+ 4H) _ 
9x2y2?BesselK [2. Vx —yly2 + (x2 —y2 2/4] v2[yl,y2] 
2m ((x1 — yl)? + (x2— y2)?)" (A+) 
y23u?,/ A BesselK [3 Vl yl)? + (2 v2? (EI vatyl.y2|_ 
4m (x1 — yl)? + (x2 y2)2)*? (A + )(A+ 2m) _ 


{ous 
T 


Again, it can be partially verified that the above expression of P satisfies the 
governing system (at all points other than the singularity) by having Mathematica® 
directly evaluate and check that 


Zx0 Pix, y] = e ae x#y. (4.10) 


The nature of the singularities of P[x,y] at x = y can be now be studied by in- 
structing Mathematica® to expand each of the four components, converted to polar 
form, in a Taylor series (see (2.16)). The result shows that 


_ fottpr ott 
Plrl= a fe | 


The comments made in Remark 4.5 also apply to the matrix P. 
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4.8 Example. By taking once more (x1, x2) = (0.5, 0.5), we can illustrate graphi- 
cally how each element of P[{0.5, 0.5}, {yl, y2}] behaves when the singular point 
is in S. The graphs in Fig. 4.3 reveal the existence of a strong singularity at 
(yl, y2) = (x1, x2) = (0.5, 0.5). Such singularities are especially troublesome in 
numerical integration, but, fortunately, this is not an issue in our case when x is off 
the boundary. 


Fig. 4.3 The components of P[x,y] x,y € S. 


The presence of a strong singularity on the boundary, however, gives rise to dif- 
ficulties. 


4.9 Example. Suppose that a segment of the line yl = y2 which contains the 
point (x1, x2) = (0.5, 0.5) is part of the boundary. The graphs in Fig. 4.4 show 
the behavior of the four components of P[x,y] in the neighborhood of the point 
(yl, y2) = (0.5, 0.5). These 3-D graphs are displayed to the left and above the 
boundary curve and truncated on the opposite side. In this case, the unit normal 


vector is v = (v1(y), v2(y)) = (1/V2, -1/v2). 
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1 


Fig. 4.4 The singularity in the components of P[x,y], x, y € OS. 


Direct computation indicates that the behavior of the components P;, g [x,y], con- 
verted to polar coordinates, near x = y on the boundary line yl = y2 is given by 


Pir] = i oll/ ) | 
O[1/r] Ofr] 
4.10 Remark. Although P, [x,y] and P,,2[x,y] contain factors that suggest the pres- 
ence of a strong singularity at x = y on 0S, analytically this is not the case, since 
the combinations of these factors with the components of the normal derivative turn 
the leading diagonal elements of P|x,y] into bounded functions when both x and y 
are on the boundary. In fact, as already mentioned above, these functions tend to 
0 at x = y with order O[r]. However, numerical techniques never calculate P|x,y] 
exactly on 0S (see Remark 4.5), and the graphs of P; 1 |[x,y] and P22[x,y] shown in 
Fig. 4.4 indicate that if y is slightly off the boundary in the vicinity of (0.5,0.5), it 
is close to a strong singularity. This, along with the strongly singular behavior of 
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the other two components, creates significant problems for numerical integration. 
The difficulty is made worse by the fact that the internal definitions of D and P are 
much more complicated than what their graphs suggest. These issues are discussed 
in much more detail below. 


4.5 The Singularities of D and P 


4.11 Example. We re-examine the behavior of D[x,y] near its singularity by con- 
sidering the earlier setup, where the singular point is (x1,x2) = (0.5,0.5) and the 
boundary curve includes a segment of the line yl = y2 containing this point. The 
values assigned for this illustration are 


A>1, w>2, k733. 


Figure 4.5 shows the graphs of the Dap [x,y] = Dap [{0.5, 0.5}, {y1, y2}] along the 
boundary segment 


yl=y2, O<yl<l. 


We notice the logarithmic singularity in D; ; [x,y] and D [x,y] and the fact that the 
other two components do not have a jump discontinuity at (yl, y2) = (0.5, 0.5). 


0.2 


. 0.02 
0 0 
0.5 1 0 0.5 1 
0.02 0.2 
0 0 
0 0.5 1 0 0.5 1 


Fig. 4.5 The singularities of D[x,y], x, y € OS. 


Since the algebraic structure of D is very complicated, it is desirable to make a 
number of simplifying assumptions so that a more detailed examination becomes 
possible near the singularity. Thus, we consider the special case where the two- 
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dimensional domain of D is rendered one-dimensional by considering D[x,y] to 
be defined on the boundary segment yl = y2 with the singularity at the point 
(x1, x2) = (0.5, 0.5) on it; that is, we make the replacements 


(x1, x2) + (0.5, 0.5), 
y2— yl, 


Vx yl)? +(x2—y2)?>r. 


This changes D[x,y] to Dir], and we can now perform a single-variable analysis 
on the new version. First, consider the dominant terms in the Taylor series expansion 
of D[r]: 


3 
5 


3 
—7EulerGamma — 2Log | w] 5Log | | TLog|r] 


407 


3 
—7EulerGamma — 2Log Se] 5Log 


a TLog{r] 


3 ik 
aon 1 ON" 40n 


+ O[r]! 

For small values of r, we have the asymptotic behavior (4.7), which confirms that 
the leading diagonal terms have a logarithmic singularity and that the other two are 
bounded as r — 0. 


It is important to bear in mind that the above results have been calculated only for 
a specific case; hence, their simplicity disguises the underlying difficulties encoun- 
tered in the construction of D. For example, D[r] contains modified Bessel functions, 
which have a singularity at r = 0: 


BesselK(0, r] > Log[r], 
BesselK{1,r] > 1/r, 
BesselK[2,r] > 1/r’. 


These functions are combined with other algebraic expressions in such a way that 
all but the logarithmic singularities cancel out in the leading diagonal entries as 
r — 0, and all the singularities in the other two components are eliminated. However, 
when D is evaluated numerically, the expressions involving D must be evaluated 
individual part by individual part before they are amalgamated into a final result. 
A detailed examination of Dap indicates that, for small values of r, the internal 
structure of D has the asymptotic behavior 


O[rJOft/r] Of JOl1/r) — Ofr2)0f1/77 
Ho Plot rio | of /r4. 


(4.11) 


Daplr] = O[Loglr]] 4 
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4.12 Remark. Formula (4.11) applies to the specific algebraic representation used 
in our code and may change if another such representation is chosen. 


This shows that, individually, the various parts of D have a much stronger sin- 
gular behavior than what we expect when they are assembled together after internal 
cancelation, which makes the numerical calculation considerably more inaccurate 
than would otherwise be the case. It is a well-understood concept in numerical anal- 
ysis that the difference of nearly equal large expressions generates a significant loss 
of accuracy when floating point arithmetic is used. The problem is exacerbated by 
the higher-than-expected order of the singularities in the individual expressions be- 
ing evaluated. The worst possible asymptotic scenario when these components are 
calculated before internal cancelation takes place is 


_ (on/r) one 
ae Ge HA : 


We can convince ourselves of the loss of computational accuracy in the neigh- 
borhood of r = 0 if we divide the components of D[r] by the first term in their 
corresponding Taylor series. Then we should have 


Di, [r] D,2[r] 
— first term first term 11 
Limit ,r>o0} = li}: 
D2 [r] D22 [r] 
first term first term 


Graphing the four components for 0.0 < r < 0.00001 shows the numerical insta- 
bility of the computation near r = 0 (see Fig. 4.6). 


1.001 1.001 + 
1 1 
0.999 0.999 | 
0 0.00001 0 0.00001 
1.001 1.001 + 
1 1 
0.999 0.999 | 
0 0.00001 0 0.00001 


Fig. 4.6 Numerical instability in the computation of D[r] as r > 0. 
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The graphs in Fig. 4.6 indicate a significant loss of accuracy for small values of 
r. In this case, the loss amounts to more than 10 digits of relative accuracy, which 
has a severe adverse effect on the numerical evaluation of integrals. In the bound- 
ary element method, we expect the largest contribution to the value of our definite 
integrals to occur in the vicinity of the singular point. Unfortunately, this is exactly 
where the greatest numerical inaccuracy also occurs. 

At a first glance, a relative error of about 0.1% should be easy to handle with 
a small increase in computational working precision. But this is not necessarily 
the case. The matrix D is used in the boundary element method as the kernel of a 
definite integral along the boundary (see (4.8)). Typically, numerical quadrature is 
of the form 


DY wif x), 
where the w; are the weights corresponding to the quadrature method and the x; are 
the quadrature points. In cases involving D, the numbers f(x;) become extremely 
large for x; near the singularity. Therefore, although the relative error in the f(x;) is 
small, the values of f(x;) when the x; are in the neighborhood of the singular point 
can significantly contaminate the result of the numerical quadrature. 


Numerical integration issues raised by singularities are discussed at length in 
Sect. 4.6. 


4.13 Remark. Numerical inaccuracy can be partially mitigated by 


(i) increasing the internal working precision of Mathematica® ; 
(ii) simplifying the expression of D to force as much internal cancelation as pos- 
sible; 
(iii) eliminating the higher-order modified Bessel functions, which contain the 
higher-order singularities. 


The default floating point working precision is machine-dependent, with a com- 
mon value of 16 digits. The graphs in Fig. 4.6 have been generated with this value. 
If the working precision is increased by 10 digits—that is, from 16 to 26—then 
most of the numerical computational inaccuracy disappears. The graphs for D[r] in 
Fig. 4.7, produced by the same algorithm but with this increased working precision, 
confirm the expected improvement. 

Simplification in Mathematica® can be performed by means of the functions 
Simplify and FullSimplify. They both operate by applying appropriate 
transformation rules to their argument and returning the simplest possible expres- 
sion for it. One problem here is that the number of possible simplification paths that 
must be explored grows tremendously fast with the size of the input expression. As 
a consequence, the process can take an excessive amount of time to reach comple- 
tion. When large expressions are encountered, such as one of the components of 
D, it is usually necessary to segment the problem into smaller pieces and simplify 
these preliminary results before simplifying the full expression. Simplify and 
FullSimplify often need additional information to complete the process. This 
is supplied by inserting suitable assumptions in the program. 
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Fig. 4.7 Improved numerical stability near the singularities of D[r] as r > 0. 


4.14 Example. The simple algebraic expression VP? that appears repeatedly within 
D{r] simplifies differently depending on the assumption on the nature of r. Specifi- 
cally, 


Vr simplifies to r if we assume that r € Reals and r > 0; 
Vr simplifies to Abs|[r] if we assume that r € Reals; 


VP simplifies to Vr? if we assume that r € Complex, which is the default. 


4.15 Remark. The function Simplify is restricted to the most commonly used 
simplification procedures, usually consisting of algebraic, exponential, and trigono- 
metric transformation rules. The expressions of D and P contain additional func- 
tions, including the modified Besse1K functions. Identities related to the Besse1K 
are not included in the structure of Simplify, so this function will not achieve the 
desired simplification. In such cases we need to use the function Full Simplify, 
which contains a much larger set of transformations but is correspondingly much 
slower in its evaluation. 


4.16 Example. The function Simplify is inadequate for 
rBesselK(0, r] — 2BesselK[1, r] + rBesselK{2, r], 

whereas Full Simplify produces the result 
2rBesselK(0, r]. 


The original expression of D is very large, with a LeafCount of 2,943. How- 
ever, after extensive simplification and forcing the elimination of the higher-order 
modified Bessel function BesselK[2,r], it can be reduced to a much smaller one, 
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with a LeafCount of only 782. Also, the internal structure of D at the singularity 
now has the more computationally friendly form 


rjo[r! 
OfLog|r|) + O[Log|r]] + O[r]O[r"] 
ae -1 O[r|O[r-} 
OlLog(r]] + Olr}Ol-'] OlLegtr]] + Gra — 


As mentioned in Remark 4.12, this expression corresponds to the algebraic form 
chosen for our code and, as such, is not unique. 

Next, we consider the algebraic structure of P, which is even more complicated 
than that of D and was examined in the previous section. 


4.17 Example. We go back to the earlier special case (see Example 4.11), where the 
singularity is at (x1, x2) = (0.5, 0.5) and the boundary curve includes a segment of 
the line yl = y2 that contains the singular point. The other parameter values are 
unchanged, and the unit normal vector is (v1, v2) = (1/v2, — 1/V2). 

The graphs in Fig. 4.8 illustrate the behavior of the components 


Pa.pl%sy] = Pa,pl{0-5, 0.5}, {yl, y2}] 
along the boundary segment 
yl=y2, O<yl<l. 


We notice the strong singularity of P;,2[x,y] and P),1[x,y], and the absence along 
the boundary path of any singularity or discontinuity at (yl, y2) = (0.5, 0.5) in the 
other two components. 
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Fig. 4.8 The singularities of P[x,y], x, y € 0S. 
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The analysis of P is in many respects similar to that of D, but there are also 
some important differences between the two. After all the internal cancelations, the 
theoretical order of the singularities in P when y lies on the boundary is 


r= (Oe on |) 


Close examination of P reveals that it contains the modified Bessel function 
BesselK(3,7], which is of order r~> for small values of r. As a consequence, P is 
constructed from parts with worse singularities than those in D. A thorough check 
of their orders before internal cancelation indicates that, as r + 0, 


O|r|O|Log|r O|r?]O[Log|r O[r\o[r—! O[r\ofr-! 
Adie [ — ]] te [ a g[r] : K : ] % a ~ | 
" O[r]O[r-|] O[rjO[r-? _O P|o[r-?| a O[r?]O[r->] 
Or) OP} Ofr4] Or] ’ 
ne O[r]O[Log[r]]_, Ofr’]O[Logir]]_, Ofrlo[r—"] ‘e O[r7]O[r-|] 
= Or] Ol] "~~ Olr] Or] 
O[r]O[r-7] O[r?]O[r-? O rjo[r-}| 
+~or) 1 Of] Or)’ 
_ Ofr|jO[Log{r]] | Of ]O[Log[r]] | Ofjo[r~"] , O[r’Jo[r—'] 
Pall = ~~ ony tal Or) | OF) 
O[rjO[r-7]__ Of |O[r-? _O r|o[r->| 
+~or) 7 of) * ory” 
[rJO[Logir]] , Ofr-]O[Log[r]] _ Ofrjofr"] _, ofjo[r—") 
Pall= “op, +o], on OF] 
O[r7]O[r-|] O[rjO[r-* of P|o[r-?| O[r-]O Fr | 
*~or) * or) ° of) ° ory 


The statement made in Remark 4.12 applies here as well. 

As in the case of D, we see that the individual parts of P have a much stronger 
singular behavior than that of the final result obtained by combining them. The worst 
possible asymptotic scenario here for the numerical singular behavior of P[r] is 


o[r“} ofr“) 
ofr] of 4]) 
To show graphically the loss of computational accuracy for P[r] as r + 0, we note 
that P; ;[r] and P)2[r] have a zero of order O[r] and are, therefore, bounded, whereas 
the other two components, which exhibit singular behavior, are unbounded. The 


latter can be divided by the singular portion of the first term in their Taylor series, 
namely 
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—(1/(5V22))r1 = Of!) for Pals), 


to generate bounded results. The expected outcome for all the components of P is 
that they should tend to either 0 or 1 as r — 0; that is, 


Pi2\r 
Palr| arial re 
Limit rstterm | +0] = 
Py [r] P. [ 1 0 
first term seca 


Figure 4.9 shows the numerical instability of the computation of the four com- 
ponents for small values of r. 


1.001 
1 
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0 0.00001 0 0.00001 
1.001 10 
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0 0.00001 0 0.00001 


Fig. 4.9 Numerical instability for P[r] as r > 0. 


The graphs in Fig. 4.9, generated for the same range of r that we used for D[r], 
indicate a much greater loss of accuracy than in the case of D[r]. This is especially 
relevant for the comparison with D, ;[r] and D2 2[r], which have a magnitude of at 
most 0.006 but produce errors with magnitudes significantly larger than 12. Such an 
extreme relative error may seem surprising at first, until we realize that D, ;[0] and 
D2.2[0] are close to an off-boundary singularity in the interior of the domain, with 
theoretical order O[r~'] and numerical order O[r~‘]. It is therefore obvious that we 
need to develop special strategies for handling numerical difficulties of this nature 
successfully. 


The procedures for reducing the numerical inaccuracy mentioned in Remark 4.13 
apply in equal measure to P. For example, Fig. 4.9 has been generated with 16-digit 
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working precision. Increasing the working precision by 10 digits causes most of the 
computational inaccuracy to disappear, as can be seen from Fig. 4.10, which has 
been produced with the same algorithm but with 26-digit working precision. 


10 1.001 | 
0 | ete 
-10 0.999 
0 0.00001 0 0.00001 
1.001 10 
ee 0 
0.999 -10 
0 0.00001 0 0.00001 


Fig. 4.10 Improved numerical stability for P[r] as r > 0. 


The matrix P can be simplified considerably. After that, and after the func- 
tions BesselK{2,r] and BesselK{3,r] are forcibly removed, we can reduce the 
LeafCount for P from 29,183 to 3,477. This shrinks the size of P by about a 
factor of ten. Also, the internal structure of the singularity in P now has the more 
computationally friendly asymptotic form 


_ OfrOfLogir]] | Of] OfLogir]] | ofjofr'} , oflor"] 
Puall=—“om, + orl + on OF 
O[r’| O[r-"} 
+ Or,” 
_ Of|OlLog(r]]_ Of} OLegir]] | oof!) _ e[r*] ofr") 
Pali=—“oy + or, + on * oF, 
py fo) = QeiolLestell , OL] olLoetr| | Oblol'} , OfF] or 
ae OP] O[r4] Or) Of]? 
Bah = OlrlO[Log{r] __ O[r'] OfLogir]] | Ofrjofr"'] _, Ofr|olr—") 
ee OF ofr Or] OFF] 
o[r] O[r-'} 
*~ OP 


However, the order of the computational singularity is still higher than the expected 
theoretical order of O[r~']. 
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This analysis is used later in the chapter to determine how to handle numerical 
integration issues that arise in the boundary element method. These issues include 


(i) the ability to handle logarithmic singularities; 
(ii) the ability to handle Cauchy principal value integration; 


(iii) the ability to increase computational accuracy to compensate for the inaccu- 
racy introduced by the behavior of D and P near their singularities. 


4.6 Numerical Integration 


Anyone who spends a lot of time working with singular integrals eventually spends 
a lot of time struggling with numerical integration techniques. In general, the algo- 
rithms coded in the NIntegrate function are near optimal for generic purposes, 
and any attempt to create one’s own algorithms usually results in disappointment. 
However, we have to do exactly that since in solving our problems we need to handle 
both weak and strong singularities. 

The first consideration is the degree of desired accuracy in the computed result. 
This is ensured by specifying the Mathematica® parameters AccuracyGoal, 
PrecisionGoal, and WorkingPrecision. Ata first glance, it might appear 
that using machine precision (as mentioned earlier, normally about 16 digits) would 
be sufficient to meet our needs. But the definite integrals of functions with singular- 
ities present special problems. Below, we show in greater detail that there is often a 
significant loss of accuracy in the computation of integrals of such functions. 


4.18 Example. The standard strategy for computing the Cauchy principal value of 
the integral over an interval (a,b) of a function f with a strong singularity at a point 
c, a<c <b, is to use the decomposition 


b c-é € b 
[r@a= fl rears [(Fe+8)+Fe-E)ds+ f far, 4.12 
a a 0 


c+E 


dictated by the necessity to integrate f on both sides of the singularity over a sym- 
metric interval of length €. The two functions f(c+ 6) and f(c—&) have nearly 
equal extremely large absolute values but are of opposite signs on 0 < & < €, which 
leads to a significant loss in precision since most of the significant digits for each 
function cancel out during the subtraction. The remaining digits are of much lower 
accuracy. The closer we are to the singularity, the more significant digits we lose. To 
compensate for this effect, an increase in computational accuracy along with helpful 
transformations must be deliberately built into the method. 


The parameters that control accuracy and precision are 


PrecisionGoal — p, AccuracyGoal — a, 
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and they tell Mathematica® to attempt to make the numerical error in a result x to 
be less than 10~¢ + |x|10~?. This would appear to be very straightforward, except 
that special problems occur in certain circumstances. 


4.19 Example. Many of the spline basis functions used to approximate the solu- 
tion are numerically zero on certain subintervals of the boundary. In these cases, it 
is clear that the PrecisionGoal for the piece of the integral over that specific 
subinterval can never be met because x is 0. A lot of computational time can be 
consumed by numerically integrating 


Tit 
i 0.0dt, 


tj 


which produces the error message “Integral and error estimates are 0 on all integra- 
tion subregions. Try increasing the value of the MinRecursion option. If value 
of integral may be 0, specify a finite value for the AccuracyGoal option.” 


Using a value for AccuracyGoal would appear to be a good solution, but it 
may be unclear ahead of time exactly what an acceptable magnitude of the definite 
integrals being evaluated might be. This can cause us to set an unacceptable error 
tolerance, either much too large or unrealistically small. The best way to handle 
this is to preprocess the definite integral symbolically and set its value to 0 on the 
relevant subinterval before it is sent to the numerical integration algorithm. 

Ill-conditioning is another issue encountered when using numerical methods. 
Frequently, the default machine accuracy of 16 digits is insufficient to overcome 
it, and we are forced to increase the WorkingPrecision option to get a satis- 
factory resolution. This, however, presents a new problem. As all Mathematica® 
programmers know, a floating point expression will automatically assume the pre- 
cision of the least precise subexpression it contains. Consequently, entering, say, 
1/2 as 0.5 eliminates any possibility of computing the result with greater accuracy 
than is available with the default machine precision. Since it may be unclear at the 
beginning how much additional precision is required to achieve our objectives, it is 
best for us to start with exact arithmetic for all subexpressions and then, if necessary, 
increase the working precision beyond the default machine precision when the com- 
putationally sensitive portion of the program is reached. Once the computationally 
intensive part has been completed, we can adjust the working precision downward 
and thus save computational time. Exact arithmetic is easy to code. For example, 
the machine accurate floating-point number 3.1415 would be entered as the infinite 
precision rational number 31415/10000. 

Another issue is overcoming the loss of accuracy in computing integrals in the 
vicinity of a singular point. This usually requires increasing the value assigned to 
the WorkingPrecision option even more. Normally, an increased working pre- 
cision is required only for the evaluation of the integral over the subinterval con- 
taining the singularity; maintaining it outside that region would significantly slow 
down the computation of all the other results. We take the view that determining 


56 4 Software Development 


the necessary additional increase in working precision is a matter of experience, as 
well as trial and error. Often, an additional 30 digits of working precision accuracy 
are required to achieve the PrecisionGoal specified. This could result in a total 
working precision exceeding 40 digits. The code used in the creation of the material 
for this book has been made adaptive in the sense that if the PrecisionGoal has 
not been met, then the current working precision is automatically increased, and the 
result is recomputed until convergence occurs with the desired accuracy. 

In summary, our objective is to take the interval (a,b) and divide it into subin- 
tervals using the procedure described above, so that f(t) is analytic in each of the 
subintervals (a;,b;); that is, 


b noIntervals 
[toa > [far 
a j=l gg 

uy: 


Next, the function NIntegrate numerically evaluates the corresponding sub- 
integrals by means of adaptive techniques, until the preset precision goal is met for 
each of them; specifically, 


Hf 


b bj 
noIntervals noIntervals order 
froa= y [roae YY waste. 
i j=l aj j=l i=1 


Special attention must be paid to the subinterval containing the singularity, which 
needs to be handled on its own. Consider, for example, a Cauchy principal value 
integral with the singularity in the integrand f at t =c, aj <c < bj. By (4.12), 


bj c-é € bj 
[roa= f soars fie+8)+Fle-8))db + f flat. 
aj aj 0 cC+E 


The value of € is chosen to be a fraction of the distance from the singular point 
to the nearest point on the Exclusions list compiled for NIntegrate, which 
comprises all the singularities, the corners that the boundary might possess, and the 
knots used in any spline approximation. 

The two integrals 


ce bj 
[ fear. [fea 


are evaluated by NIntegrate in the usual way. The singular integral 


E 


[ule+8)+Fle-S) a8 


0 


is now evaluated separately by NIntegrate using increased working precision 
combined with singularity-handling adaptive numerical quadrature techniques. In 
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our experience, the double exponential quadrature method yields the fastest and 
most accurate result for the type of singularities encountered in the boundary ele- 
ment technique. This method transforms the original integrand in a way that sup- 
presses the singularity and subsequently evaluates the integral numerically using the 
trapezoidal rule. The transformation used is specifically chosen because its deriva- 
tive decreases doubly exponentially when the variable of the integrand reaches 0 and 
é. In our case, a double exponential transformation of the form € = @(7) changes 
the integral 


0 
where 
a(g)=fle+S)+-f(e—¢), 
to 
J so(n))o' nan. (4.13) 


The integrand g(€) must be analytic in (0,¢) and may have a singularity at 0. 
The analyticity condition is satisfied because of the way we have constructed the 
Exclusions for NIntegrate. The singularity in g() must be weaker than 
O(E~') for the integral to converge. Here, the function g(E) = f(e +€)+ f(c—&), 
which cancels the strong O(€ ~') singularity, still has a logarithmic singularity of or- 
der O(Log|/&]) as € — 0. The logarithmic singularity can be compensated for by the 
use of the double exponential quadrature method because the transformed integrand 
decreases doubly exponentially; that is, 


@(6(m))6"(m) © exp(—cexp(|n]|)) as 1 > +e. 


The function (7) is analytic in (—c°, ce), The trapezoidal rule was chosen be- 
cause it is known to be optimal for this type of situation. The transformation used in 
our case, where the singularity is located at € = 0, is 


& = 4e+ 4 etanh (5 msinhn) (4.14) 
and converts € = 0 into n + — and € = € into n — »~. The trapezoidal rule applied 


to (4.13) yields 


[ s@no'nan =n Y Flo)", 


——) 


The terms in the trapezoidal rule decay doubly exponentially as i + +o; there- 
fore, the infinite sum in it is truncated when the computed terms become too small 
to make a significant contribution to the total. 
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4.20 Example. The construction of P[x,y[t]] usually guarantees that a logarithmic 
singularity still remains even after the cancelation of the strong singularity in the 
Cauchy principal value integration method. Suppose that we want to evaluate 


1 
[roelélaé. 
0 


The graph of the integrand Log|&] for 0 < € < 1 is on the left in Fig. 4.11, and the 
graph of the function (4.14) for —co < 7 < © is on the right. 


Fig. 4.11 Graphs of Log[€] and § = 5 €+ 4 €tanh (5 msinh1). 


The transformed integral is 


J so(n))o'n)an = f Loe(o(n))o'(n)an- 
The function @’(17) for € = 1, which contains the double exponential decay, is 
¢'(n) = } aCosh{n]Sech[} mSinh{n]]’, 


graphed on the left in Fig. 4.12. We remark that this function decays very quickly, 
as seen from the values 


o'(+3)=6.8x 1078, '(44) =5.0x 10-6. 


Fig. 4.12 The double exponential decay transformation. 
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The graph on the right in Fig. 4.12 represents the integrand 


g(o(n))¢"(n) 
= (Log[4 +! Tanh[4 Sinh{n]]]) ({ tCosh{n|Sech[} mSinh[n]]°). 


We notice that 
g(o(3))o’(3) = —2; x 1074 «6.8 x 10719, 


and that the logarithmic singularity has been completely overpowered. 


Example 4.20 illustrates three important features of the double exponential 
quadrature method: 


(i) The domain of numerical integration is finite, not infinite, and is usually quite 
narrow. 

(ii) Additional digits of accuracy can be achieved with a very small increase in 
the size of the domain of numerical integration. 


(iii) The singularity in the integrand has been removed. 


The first feature implies that the number of functional evaluations required in the 
trapezoidal rule can be reduced, which increases the computational speed of the 
method. The second one means that adaptive quadrature requires very little addi- 
tional work to increase its level of accuracy. Finally, the third feature points to the 
fact that, the singularity having been removed, the trapezoidal rule will not involve 
operations with numbers of significantly different magnitudes. Combining large and 
small numbers in floating-point arithmetic usually results in a slowing of the rate of 
convergence as well as a deterioration in accuracy. 

The double exponential quadrature method can be much faster than other quadra- 
ture techniques. 


4.21 Example. The regular (untransformed) trapezoidal rule has a great deal of dif- 
ficulty with the singularity in an integral such as 


1 
[bosiglaé. (4.15) 
0 


We can compare the two methods on the basis of the number of digits of accuracy 
required in the answer. The bottom set of points on the left in Fig. 4.13 shows that to 
achieve | to 5 digits of relative accuracy, the double exponential quadrature method 
needs only about 12 to 30 evaluations of Log[&]. The top set of points shows that 
the untransformed trapezoidal rule requires over 1,000,000 evaluations of Log[é] 
to ensure 5 digits of accuracy. Clearly, the time consumed by the untransformed 
trapezoidal rule in a case like this is unacceptably excessive. 


4.22 Example. We can also compare the double exponential quadrature method to 
the global adaptive technique, which is the default procedure used by NIntegrate. 
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The latter is a higher-order Gaussian quadrature method with adaptive strategies and 
singularity-handling capabilities, and it does a reasonable job in evaluating the inte- 
gral (4.15). 
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Fig. 4.13, Double exponential quadrature vs. global adaptive method. 


The graph on the right in Fig. 4.13 shows the number of evaluations of Log[&] 
required to achieve 50 digits of relative accuracy for each quadrature method. We 
see that the double exponential quadrature method (the bottom set of points) does 
this with about 300 evaluations as opposed to over 2,000 evaluations needed by the 
global adaptive method (the top set of points). 


4.23 Remark. It is useful to know where and how often the integrand Log{&] in 
(4.15) is evaluated to perform the numerical quadrature. The graph on the left in 
Fig. 4.14 shows the exact sequence of chosen values of €, 0 < € < 1. The desired 
relative accuracy for our answer was 30 digits. We notice that more evaluations 
occur at the endpoints of the domain, where we anticipate that a singularity might 
exist. 
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Fig. 4.14 Double exponential quadrature. 
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The graph on the right shows more clearly the integration scheme in the trans- 
formed variable 1, —co < 1 < oo. First, we see the limited range of 7 values actually 
used. Also, we notice the implementation of the trapezoidal rule in the transformed 
variable. This is recursive, with the number of evaluation points doubled every time 
the accuracy is insufficient, until convergence occurs. 


Next, we investigate how the double exponential quadrature method changes the 
behavior of D in the vicinity of its singularities. To this end, we use again the form 
D[r], 0<r< 1, derived for D in Example 4.11. The loss of relative accuracy in the 
computation of the four components of D[r] can be seen from the graphs for 


0.0000001 < r < 0.00001, 


shown in Fig. 4.15. These graphs indicate numerical instability as r — 0. They were 
generated with a working precision of 16 digits, and the range of 1 + 0.001 repre- 
sents a relative error of +0.1%. 


1.001 1.001 + 
1 1 
0.999 0.999 + 
0 0.00001 0 0.00001 
1.001 1.001 + 
1 1 
0.999 0.999 | 
0 0.00001 0 0.00001 


Fig. 4.15 Relative computational error for the untransformed D[r] as r > 0. 


The double exponential quadrature method transforms the integral 
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where the integrand, when generated by the double exponential transformation 
(4.14), is 


D(o(n))¢'(n) 
= (D[4 +4 Tanh[+ rSinh[n]]]) (£ tCosh{n]Sech[4 rSinh{n]]”). 


The graphs in Fig. 4.16 show the behavior of the transformed matrix D(@(n))’(n). 


0.06 ; 


0.03 | 


—2 -l 0 1 2 


Fig. 4.16 The transformed matrix D(g(1))@'(7). 


These graphs have been drawn for 
—234< 1 <-2ol, 

which corresponds to 
0.0000001 < r < 0.00001. 


We recall that D[r] is not bounded as r — 0, but that its relative size is bounded. The 
transformed expression D(@(7))@’(1) is bounded everywhere. 

Dividing each of the four components by their maximum value, we can plot the 
relative error for the transformed expressions. The graphs in Fig. 4.17 show the 
relative error in the computation of D(@(7))‘(1) over the equivalent domain in 7; 
the range (0.0, 0.001) represents a relative error of 40.1% of the maximum value 
for each component. These graphs, generated with a working precision of 16 digits, 
do not exhibit any of the numerical instability seen earlier in the untransformed D[r] 
for small values of r. 
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Fig. 4.17 The relative computational error in the transformed D(g(n))o'(n). 


In summary, the double exponential quadrature method represents an ideal choice 
for performing the numerical integration required by the boundary element method. 


4.7 Coding the Single-Layer and Double-Layer Potentials 


With the necessary symbolic code now in place, we can formulate the single-layer 
and double-layer potentials (see Definition 2.3) as 


-f Pix, y]o wly] (4.16) 


where @ = {1 [y], @2[y]} and w = {wi|y], Wo[y]} are two-component density func- 
tions. When x is on the boundary, D[x,y] and P[x,y] have weak and strong singular- 
ities, which the code designed to evaluate the boundary integrals containing these 
matrices as kernels needs to take into account. When x is not on the boundary, the 
kernels have no singularities, and normal integration along the boundary can be 
performed to evaluate the corresponding boundary integrals. With this coding, the 
Somigliana formula (2.18) with x € OS is written as 


¢ Deso(T)bld= f Plxyloublar+ hubs, (4.17) 
TWeak{x] TCPVv ix] 


and with x € S, 
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~o, 
S 
x, 
(e} 
S 
= 
ra 
lI 


f Plxy]oub dI, +ulx]. (4.18) 
T T 


4.8 Coding the Boundary Integral Methods 


After the single-layer and double-layer potentials have been coded, it is straightfor- 
ward to do the same for each of the boundary integral methods developed in Chap. 
3. For brevity, we illustrate only the case of Dirichlet boundary data. 


4.24 Example. The direct method for the Dirichlet problem (D*) defined in Sect. 
3.2 leads to the boundary integral equation (F; ); that is, 


Vo(~) = (Wo+ 51) F, 


where # is the prescribed Dirichlet data function on 0S. Using (4.17), we code this 
equation as 


$ Dix, y] o gy] - ¢ Pix,y]o Aly dTy yt 5 Ala, (4.19) 
Fweakis| Topvpy 


where @ is the unknown density representing Tu. Once this density is determined, 
we can find the solution in S by means of the representation formula (4.18): 


= p Dix.yo@bldly— $PxyloPhldh, xes. (4.20) 
EF Tr 


4.25 Example. The classical indirect method for the Dirichlet problem in S, de- 
scribed in Sect. 3.1, leads to the boundary integral equation (ae ); that is, 


(Wo- 31) 0=F, 


where # is the prescribed boundary data function on 0S. This equation is coded in 
Mathematica® as 


TCpVv ix] 


where @[y] is the unknown density. Again, once the density is determined, we can 
find the solution in S by means of formula (4.16). 


The coding shown above omits a lot of detail. Basically, code is developed to 
convert the boundary integrals into parametric line integrals. The latter are then 
evaluated numerically, special attention being paid to the boundary singularities of 
their kernels. 
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4.9 Outline of the Boundary Element Method 


The boundary element method uses certain functions (elements) to approximate the 
unknown densities of the boundary integral equations, and consists of several steps. 


Step 1. We write the symbolic quantities D and P corresponding to the homoge- 
neous problem 


Zyou=0. 


Step 2. We construct the single-layer and double-layer potentials 


(Ve)lx| = a Dix.y]o@bld 


Topy [x 


Step 3. We code the equation for the specific boundary integral method we want to 
apply. Thus, equation (ft ) in the direct method for the Dirichlet problem takes the 
form shown in Example 4.24, namely 


$ Dix, y] o gy] - 4 Plx,ylo Aly| dT, yt5 Al, 


TF Weak Tepvis] 
where @ is the unknown density and # is prescribed on 0S. 


Step 4. We define a parameterization y|t] of the boundary 0S. This is done by split- 
ting the boundary into segments (0S); corresponding to the Dirichlet, Neumann, or 
Robin conditions prescribed on them, respectively, or to any corners the boundary 
might have. Corners introduce discontinuities in the derivative of the parameteri- 
zation and, therefore, create potential discontinuities in one or several of y’ |r], v[t], 
Dix, y[t]}, and Plx,y{¢] 


Step 5. We use the parameterization of the boundary to change the boundary integral 
equation into a line integral of the form 


D{x, y[t}] 0 p[y[t]]dsdt/r] de 
TWeak{x] 
$ Px, yl] o Ply le]dsdt{e] dt +4 Ply, 
TOpvpy 


Step 6. We construct approximate solutions @ by using the parameterization y/f] 
and choosing the knots ¢;; and the spline degree required to make up a B-spline 
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basis {b;,;}. The location of the knots is dictated by the boundary segments (0S);. 
Additional knots within each (0S); are also selected in order to achieve the desired 
approximating properties. The approximate solutions are represented as linear com- 
binations of the b;,; with unknown coefficients Cg ;,;; that is, 


~ (dedi ji,j 0 
O[y[]] = &S a , 


Step 7. We substitute the spline elements in the boundary integral equation. For 
example, equation (Ds ) in the direct method for the Dirichlet problem becomes 


fi rows & ee asa 
J Sab ITU 


TWeak( 
~ $ Plx,y[t]]o Ply|eljdsat[e] dr + 1 Ply). 
TCPVvix| 


This converts our boundary integral method, which is exact, to a boundary element 
method, which is approximate. 


Step 8. Because Mathematica® is a list-based language, the bookkeeping does not 
become an issue. The list containing the B-spline basis {b; ;} for each of the two 
components of @ can be combined into a single list twice as long; more precisely, 


b;,;{t] 0 
hay 0 seri bi. i(t] vee 


is converted to 


LelDpg[t|| a 0) 0 0 
0 0 0 ... Blt]... : 


The 2 x 1 matrix equation of the boundary element method at the point x can now 
be written as 


ae Cli,j 
$ D{x,y[t]] (; uae . 0 me ") o| ... | dsdt(t]dr 
TWeak(x} ee : C23; 
~ f Pix, y{t]] 0 Plylél|dsat(e] dr + 4 Aly. 


TCpvix| 


Using the linearity of the integral, we bring this equation to the form 
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C1ij 


( f Dix yfelJo i oe . . oo °) dsdt(t] ar) ° 


TWeak([x] C24, 


~ ¢ Ple.ylllo Ppleldsadl dr +4 Ply. 


TCPVvix| 


Step 9. We select a method for computing the unknown coefficients cg,;,;. Although 
there are several methods to choose from, here we use only the collocation tech- 
nique. We start by selecting the collocation points x,, k = 1,...,m, on dS and map 
the 2 x | equation of the boundary element method onto these points to form a linear 
system of 2m equations. The x, must be chosen so as not to include any corners, as 
well as to ensure that the resulting linear system has a solution. This system is of 
the form 


ae 0 0 0 
D[x,,y[t}] 0 iM dsdt[rJdr }° |: 
0 0 bj. ;\t 
R a C2.) 
Weak [xj] . 
= $ Plxe,yltl]o Pplelldsat(e] a+} P [x] 

Topv [x] 


The x; also include the location of the singularities in D[x, y[t]] or P|x, y[t]]. These 
points must be explicitly passed to the numeric integration algorithms in order to 
ensure accurate evaluation. 


Step 10. The individual components of both the coefficient matrix and the right- 
hand side vector above are sent for numerical integration. The integration algorithms 
must accurately handle the singularities of the integrand. This process returns a 
system with numerical coefficients. 


Step 11. The system constructed in the preceding step is either ‘square’ or overde- 
termined, and is solved for the unknown coefficients cg, ;j, ;. 


Step 12. The approximate solution for the unknown density is now represented in 
terms of the spline basis elements as 
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Oly[t]] = a ee 4 iis - Did 14, j5i,j[E] 
= 0 0 O ... bi ;[t] wis : Did; c24jbi,j[C] . 


Step 13. This approximation is used to determine the approximate solution wu in S. 
In the direct method for (D~), the representation formula (4.18) yields 


Cli 


ux] © (f Pix rtile & we *) Asati a) : 


r C23, 


7 f Pix, y[t]] 0 Ply[e]}dsat{e] ar. 
Fr 


The rest of this chapter provides additional details for steps 4, 6, and 8. 


4.10 Parametrization and Segmentation of the Boundary 


Before integrals over 0S can be evaluated, the boundary curve must be parameter- 
ized so that the integrals are converted to ordinary definite integrals with respect to 
a parameter ft. This is done by introducing two suitable functions 


{yl,y2} = {yl[¢], y2[]}, O<t<tMax. 


Frequently, these functions are defined piecewise, especially for a boundary with 
corners. After the parametric form for the normal {v1[¢], v2[r]} and for a [t] are cal- 


culated, the integrals over 0S are ready for conversion. For example, the Somigliana 
formula (4.17) with x on the boundary becomes 


$ Dix, ty! [¢], y2[¢]}] o (Tu) [y1[e], y2[e]] * dsdt[e] a; 
TWeak(x] 
= $ Pix, {y1[¢], y2[¢]}] oly [e], y2[e]] + dsde[e] dT; + 5 u(x]. 
Topviax| 
Numerical evaluation of integrals over 0S also requires the boundary to be parti- 


tioned into segments on which the integrand is either ‘well behaved’ or has a well- 
defined singularity. If the problem has mixed boundary conditions, the boundary 
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curve must first be divided into arcs with only one type of (Dirichlet, Neumann, or 
Robin) condition prescribed on each of them: 


as= U (dS). 


i=1,..,m 


The boundary integrals need to be evaluated numerically on each of the (0S);. Al- 
though the algorithms contained in the NIntegrate function are adaptive, they 
are based on the assumption that the integrand is nearly analytic on each segment. 
Consequently, the (0S); must be further subdivided into sub-segments determined 
by the points where the integrand is not ‘well behaved’. This happens in our prob- 
lems when D{x,y[t]] or P[x,y[t]] or their derivatives have a discontinuity for x on 
the boundary. Both u[y(t]| and D[x,y[t]] are continuous for y[f] on 0S except at the 
singularity y = x. Additionally, if the domain has corners, the derivatives of these 
functions have discontinuities at the corner locations because y’|t] is discontinuous 
there. Kernels such as (Tu)[y[¢]] and P[x, y[t]], whose structure includes the normal, 
have, therefore, a lower degree of smoothness. Hence, corners must be excluded 
from the process because the integrands are not sufficiently smooth at these points 
to meet the conditions for effective and fast numerical integration. 


4.26 Example. Let S be the upper half of the disc with center at (1,0) and radius 
1, shown on the left in Fig. 4.18. Its boundary, which consists of a half circle and a 
straight line segment, has two corners and can be parameterized by 


1 0<r<1, 
A= 
1+Cos[z(¢—-1)], 1<t<2, 
4.21 
0, 0<t<1, ¢ ) 
y2[t] = 
Sin[x(t—1)], 1<t<2 


This parameterization starts at the origin and proceeds counterclockwise around 
the boundary. Here, D[x, y[f]] is continuous for all 0 <t < 2, but its derivatives are 
discontinuous at the corners—that is, for t = 0, 1, 2. On the other hand, P|[x, y[t]] is 
discontinuous at these points, which are shown on the right in Fig. 4.18. 


OFS _ 0,@ @ 
0 1 2 0 1 


Bock 


Fig. 4.18 Left: the boundary 0S. Right: the corners and exclusions for a half disc. 
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4.27 Remark. The boundary element method requires that either u[y] or (Tu)[y] be 
replaced by approximating functions wly] or @[y| in the Somigliana formula. For 
this purpose, our approximating functions are splines, one for each of the two com- 
ponents. Splines can exhibit discontinuities, either in themselves or in their deriva- 
tives, at various knot locations, which means that the knots must also be included in 
the creation of boundary sub-segments. 


In summary, the boundary must be divided into segments determined by the na- 
ture of the conditions prescribed on it, the points of singularity of the integrands, 
any possible corners, and any knots occurring in the spline approximation. Once all 
this information is listed in the Exclusions option, the NIntegrate function 
automatically performs the segmentation of the boundary. 


4.11 Construction of a B-Spline Basis 


Suppose that the boundary 0S has been partitioned into segments in the manner 
indicated above. We are using B-splines to represent the two components of the un- 
known densities @ and y that approximate Tu and u, respectively, in the Somigliana 
formula. If we have mixed boundary conditions, then @ and yw have support only 
on the portion of the boundary associated with its specific type of condition. As 
described earlier, we expect y to be continuous on the boundary but to have dis- 
continuous derivatives at the corners. Also, we expect @ to be discontinuous at the 
corners, where the unit normal vector is undefined. 

The splines are constructed in the usual way as piecewise polynomial functions 
over a domain 0 < ft < tMax. A knot ¢; is the location between any two adjacent 
polynomial segments. At each knot, a smoothness condition can be imposed on the 
adjacent polynomials. It is not necessary to have the same smoothness condition at 
every knot. In fact, in our applications, we typically have lower smoothness con- 
ditions at the knots associated with boundary corners than at the other knots. We 
refer to the values f; corresponding to the points that separate the segments (0S); 
as primary knots. All the other knots, lying between the primary ones, are called 
secondary knots and are denoted by ¢;,;. The space between the primary knots is 
determined by the parameterization of each (0S);. The space between consecutive 
secondary knots is assumed to be equal on each individual (0S);, but not necessarily 
the same for different (0S);. 


4.28 Example. Consider a piecewise linear spline defined on the intervalO <t< 1 
with knots at 0, 1/3, 2/3, 1. In Mathematica®, specification of the knots requires 
both specification of their location and of the degree of smoothness to be imposed 
on adjacent polynomials. The function that generates the B-spline basis is 


BSplineBasis[{degree, knots},i,x], i=0,...,k. 


Here, the degree is 1, and the knots are defined with location and smoothness by the 
set 
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{0, 0, 7 3, 1, 1}. 


Figure 4.19 shows the graphs of the four piecewise linear basis functions. 


0 1/3 2/3 1 0 1/3 2/3 1 


0 1/3 2/3 1 0 1/3 2/3 1 


Fig. 4.19 Piecewise linear basis functions. 


4.29 Example. Consider piecewise cubic splines that are three times continuously 
differentiable at the secondary knots, generated with degree 3 and knot locations 
and smoothness conditions represented by the set 


{0, 0,0, 0, 5, 3, 1,1, 1, 1}. 


The graphs of the 6 basis functions are shown in Fig. 4.20. 


1 1 
0 13 2/23 1 0 13 2/3 1 
1 1 

0 1/3 223 1 0 13 2/3 1 
1 1 

0 13 223 1 0 13 2/3 1 


Fig. 4.20 The basis functions for a piecewise smooth cubic B-spline. 


4.30 Example. We now consider a piecewise cubic spline that is continuous but 
with discontinuous first-order derivatives at each knot. Its basis is generated with 
degree 3 and knot locations and smoothness conditions given by the set 
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1211222 
10;0;0, 0; ey e5es4s gags lp dy yp (4.22) 


The graphs of the corresponding 10 basis elements are shown in Fig. 4.21. 


0 13 2/3 1 0 123 2/3 1 


0 13 2/3 1 0 13 2/3 1 


Fig. 4.21 The basis functions of a piecewise continuous (with discontinuous derivatives) cubic 
B-spline. 


4.31 Remark. The bookkeeping for specifying the smoothness condition at a knot 
may appear to be a little strange, but it does have a logical explanation. A B-spline 
basis function is defined to be 0 outside the interval determined by its knots. Thus, 
in Example 4.30, each basis function is 0 for t < 0 and t > 1. The smoothness con- 
dition number m specifies the level of smoothness at the point of contact of adjacent 
polynomials and does not have to be the same at all the knots. To be specific, m = 0 
corresponds to discontinuity of the function, m = 1 corresponds to continuity of the 
function and discontinuity of its first-order derivative, m = 2 corresponds to continu- 
ity of the function and its first-order derivative and discontinuity of its second-order 
derivative, and so on. The smoothness condition number is computed from the no- 
tation for the individual knots by means of the formula 


m = (degree + 1) — (nknot repeats), 
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where ‘degree’ is the degree of the adjacent polynomials and n is the number of 
times the knot location is repeated in the knot list. 


4.32 Example. For the piecewise linear spline considered in Example 4.28, we have 
degree = 1 and n = 1; that is, each of the interior secondary knots 1/3 and 2/3 is 
listed only once. Therefore, m= 1, which means that the adjacent polynomial basis 
functions are continuous but have discontinuous first-order derivatives at these knot 
locations. 


4.33 Example. The piecewise cubic spline in Example 4.29 has degree = 3 and 
n= | at the knots 1/3 and 2/3, which yields m = 3. This means that the adjacent 
polynomials at these knots are continuous with continuous first-order and second- 
order derivatives, but discontinuous third-order derivatives. 


4.34 Example. In Example 4.30, we have degree = 3 and n = 3 at the knots 1/3 
and 2/3, so m = 1, indicating that the adjacent cubic polynomials are continuous 
but their first-order derivatives are discontinuous at the two interior knots. 


4.35 Remark. The starting and stopping knot locations require special treatment. 
As already mentioned, a spline is defined to be 0 before the first knot and after 
the last knot. In Example 4.34, each of these knots are repeated four times, which 
produces the smoothness condition number m = 0. Consequently, the first and last 
basis functions have a jump discontinuity between 0+ and 1—. The smoothness 
condition number for the second and second-last splines is m = 1, so each of them 
starts and stops with the same value (in this case, 0) but with different values for their 
first-order derivatives. Continuing this pattern, the smoothness condition number 
m steadily increases at the starting and stopping knot locations until it becomes 
numerically equal to the polynomial degree. This behavior is illustrated in Figs. 
4.19-4.21. 


We now turn our attention to the unknown densities @ and y on the boundary. 
They are approximated by 


ole] = gly[t]] Sb], O<t < tMax, 
vit] = vbte]] Sol], O<t < tMax, 


where b is a B-spline representing an approximation with primary knots ¢; deter- 
mined by all the boundary segments (0S);. Any expected smoothness in @ and y 
should be reflected by the primary knot locations t; and the smoothness at these lo- 
cations. As mentioned earlier, it is computationally unacceptable to assume a higher 
level of smoothness in an approximation than the actual level exhibited by what is 
being approximated. If the boundary has corners, then there will be discontinuities 
in the densities themselves or in their derivatives. This, too, should be reflected in the 
B-spline primary knot locations and the smoothness conditions at these locations. 
The equally-spaced secondary knots ¢;,; and their corresponding piecewise poly- 
nomials are usually chosen to achieve the approximation accuracy desired in the 
problem on each boundary segment (0S);. We recall that the density is assumed to 


74 4 Software Development 


be analytic on each (0S);. Therefore, both the degree of the piecewise polynomials 
and the smoothness at the secondary knots is dictated strictly by computational and 
accuracy considerations. 


4.36 Example. Let /: be the spacing between adjacent secondary knots ¢;,;. Then the 
expected optimal accuracy for a piecewise linear B-spline approximation is O(h”), 
and the expected optimal accuracy for a piecewise cubic B-spline approximation 
is O(h*). Whether or not this optimal accuracy can be achieved depends on the 
specifics of the selected method. 


In conclusion, the bookkeeping for constructing a B-spline from a B-spline basis 
consists of five steps. 


Step 1. The boundary OS is divided into segments (0S); with primary knots f;, 
i=1,...,m. 


Step 2. Each segment (0S); is subdivided by means of secondary equally-spaced 
knots ft, j- 


Step 3. The smoothness at each primary knot ¢; is chosen to correspond to the 
smoothness of @ or y. 


Step 4. The polynomial degree and smoothness are chosen for the secondary knots, 
thus defining the piecewise polynomial spline elements on each (0S);. 


Step 5. The basis elements 5; ; are constructed for each of the two components of 
@ or w. 


The procedure for generating the exact basis elements is somewhat complicated 
to describe and is best illustrated by a few specific cases. 


4.37 Example. Consider the half disc domain discussed in Example 4.26, with the 
boundary, corners, and exclusions shown in Fig. 4.18. The parametrization of the 
two components Wq[y[t]] of the unknown density w[y|t]], which represents u[y[f]], 
has two corners, so it produces three primary knots t; = 0, 1,2, where w[y(f]] is 
expected to be continuous but with a discontinuous first-order derivative. We intro- 
duce two secondary knots in each of the two boundary segments (0S); and (0S)2, 
and create piecewise cubic splines that are twice continuously differentiable at these 
knots. Thus, the knot locations and smoothness conditions are given by the set 


{0, 0, 0, 0, i, z, 1, 1, 1, z, 2,2, 2, parle 


Figure 4.22 shows the 11 piecewise cubic B-spline basis elements b;,; constructed 
in this way. 

In reality, the density y has periodic smoothness between the first and last basis 
functions rather than the discontinuity indicated above. Consequently, the first basis 
function needs to be modified to meet this smoothness condition, leaving only 10 of 
the b; ;, as shown in Fig. 4.23. 
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1 1 
0 18 28 1 48 53 2 0 18 28 | 48 53 2 
1 1 
0 18 28 1 48 53 2 0 18 28 | 48 53 2 
1 1 
0 18 28 1 48 53 2 0 18 238 | 48 53 2 


1 


i=) 


1/3 2/3 1 4/3 5/3 2 


1 


0 18 28 1 48 53 2 
0 18 28 1 48 58 2 


i=) 


13 23 81 48 5/73 


N 


0 1/3 2/3 1 4/3, 5/3 2 


Fig. 4.22 Piecewise continuous, with discontinuous first-order derivatives, cubic (non-periodic) 
basis functions. 


4.38 Example. With the same setup as in Example 4.37, we can also approximate 
the unknown density @, which represents Tu. Proceeding as above, we conclude that 
here the knot locations and smoothness conditions are specified by the set 


{0,0,0,0,4, 2, 1,1,1,1, 3, 2,2, 2, 2,2}. 
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1 1 
0 18 28 1 48 53 2 0 18 28 | 48 53 2 
1 1 
0 18 28 1 48 53 2 0 If 28 | 48 53 2 
1 1 
0 18 23 1 48 528 2 0 18 28 | 48 53 2 
1 1 
0 18 28 1 48 53 2 0 If 238 | 48 53 2 


0 18 28 1 48 58 2 0 18 23 1 48 53 2 


Fig. 4.23 Piecewise continuous, with discontinuous derivatives, cubic (periodic) basis functions. 


Figure 4.24 shows the 12 piecewise cubic B-spline basis elements b;,; generated 
in this case. 

The graphs of these functions, restricted to the boundary and written as b,,;[y[f]], 
are displayed in Fig. 4.25. 


After the appropriate B-spline basis functions b; ; have been constructed, we can 
define the approximate densities w and @ as 


{EE ELauvut| , O<t<tMax, 
ij J 


i i 


4.11 Construction of a B-Spline Basis 77 


1 1 


i=) 


13 2/32 1 4B 5B 


N 
i=) 


13 23 1 48 5/3 


N 


i=) 
N 
i=) 
N 


13 23 061) 48) 5B 13 23 1 48 5/3 


1 1 


0 1p 2B 1 48 5B 2 0 If 2B 1 48 53 2 
bk. Lon 
0 Ip 2B 1 48 #53 2 0 Ip 2B 1 48 53 2 


o 


1/3: 2/3 1 4/3, 5/3 


N 
i=) 


1/3 2/3 1 4/3, 5/3 


N 


i=) 


1/3 2/3 1 4/3, 5/3 2 0 1/3 2/3 1 4/3, 5/3 


N 


Fig. 4.24 Cubic B-spline basis functions for approximating @. 


= {Eesti ELeunvtl| , O<t<tMax. 


i i 


A boundary element method is then applied to solve for the unknown numerical 
coefficients Cy j,;. 
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Fig. 4.25 Cubic B-spline basis functions in 3-D. 
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4.12 Selection of the Boundary Element Method and Collocation 
Points 


The boundary integral method can be combined with the B-spline basis elements 
bj,; to produce a boundary element method. 


4.39 Example. Consider the direct method for the interior Dirichlet problem, which 
leads to equation (23 ), combined with the B-spline basis elements developed in the 
previous section; that is, 


C1ij 


( $ Dix, [tI] 0 6 Pill - bl _) dsdt[] ar) 0 
| 


TWeak(x C2,i,j 


- $ Plx,ylt]]o Ply[éljdsat{e] de + 4 Ae) 


Topvpy 


This system generates an approximation @ to the actual density @ by means of the 
formula 


wan (ADEM) — (Lider sil - 
gt] = o{t] = ee = Coe , O<t<tMax. 


Writing equation (Y;) for @ in the form 

Vo()(x) — [Wo(P)(x) +3 A(x)] =0, xe dS, 
and substituting @ in it, we define the residual 

R(x) = Vo(G)(x) — [Wo FP) (x) +5 PAQ)], xe 9S, 
which is only approximately 0; that is, 

R(x) =0, x€ds. 


We must now select a strategy for computing the unknown coefficients cg; ;. The 
usual procedure is to implement a method that makes the residual R(x) small. This 
can be accomplished in many different ways, each technique having its advantages 
and disadvantages. Our preference is known as the weighted residual method, which 
we briefly outline below. 

We start by considering a set of linearly independent functions 


{wi(x),w2(x),-..,walx)}, xe dS, (4.23) 
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that are sufficiently ‘well behaved’ to allow the evaluation of the integrals 


f R(xwile) ds, i=1,2,...,n. 
os 


The method of weighted residuals determines the cg,;,; by requiring R(x) to be or- 
thogonal (with respect to the L? inner product) to 


span {w1(x),w2(x),...,Wn(x) }. 


Thus, we want to have 


f R(owila) ds=0, i=1,2,...,n. (4.24) 
os 


These conditions generate n equations which, for the method to be successful, 
should form a consistent algebraic system. If mn equals the number of unknown co- 
efficients Cy ;,;, we may be able to find a unique solution. If the system is overdeter- 
mined, then we need to seek a least squares fit to the unknown coefficients. 

The next step is to select the functions (4.23), which immediately raises two 
issues. 


Issue 1. In the boundary element method, the residual R(x) is usually very expen- 
sive to evaluate numerically. If the w; have significant nonzero support on 0S, then 
the numerical evaluation of the equations generated by (4.24) is also very expensive. 


Issue 2. The overall accuracy of the boundary element method depends on a number 
of factors, including 


(i) the choice of boundary integral method; 
(ii) the approximating properties of the chosen B-spline basis; 
(iii) the approximating properties of the chosen functions w;; 
(iv) the choice of weighted residual method; 
(v) the accuracy of the numerical methods used in the implementation. 


Issue 2 indicates that there can be a massive number of possible variations to con- 
sider. One popular choice, known as the Galerkin method, selects the weighting 
functions w; to be the same as the approximating basis, which in our case is a B- 
spine basis. This partially addresses issue 1 because the basis has been constructed to 
have minimal nonzero support, thus simplifying the evaluation of the inner products. 
The Galerkin method with B-splines also has very good approximating properties, 
which answers issues 2(ii) and 2(ii1). 

However, issue | is the most important for our applications. Again, the expensive 
numerical evaluation of R(x) dictates the choice of the w;. Since the Dirac delta 
distribution 6(x;,x) with x; € OS minimizes the evaluation of R(x), we make this 
choice for the w;, so the weighted residual method reduces to 
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or 
R(x) =0, i=1,2,...,n. 


This procedure is known as the collocation method and is widely used in practice. 

Issue 2(iii), concerning the approximating properties of the 6(x;,x), involves both 
the number and location of the x; on 0S. It might appear that the x; can be chosen 
without much restriction, but this is not the case. The selected points must, at a 
minimum, meet the following constraints: 


(a) They must ensure that R(x;) is both uniquely defined and ‘well behaved’. 
Consequently, it is usually unacceptable to have the x; located at a boundary cor- 
ner, where the unit normal vector is undefined and the method potentially breaks 
down. There are procedures to work around this situation, but they require special 
consideration. 


(b) The choice of the x; must guarantee sufficient independent constraints to de- 
termine the cq;,;. For example, the B-spline basis is composed of functions bj; 
constructed with minimal support on the boundary. These functions have only a 
limited amount of smoothness linking adjacent polynomials. As a consequence, it is 
possible to make a selection of the x; that does not yield all the coefficients although 
the total number of these points would seem to be sufficient. 


(c) The number and location of the x; also influence the accuracy of the overall 
method. For example, their location affects the region of the boundary where R(x) 
is usually the smallest and, hence, the overall approximation accuracy. One possible 
solution to consider is selecting the x; to be the Gaussian quadrature points on each 
interval between adjacent secondary knots. Another is to select an excess number of 
collocation points x;. In practice, this choice often helps increase the accuracy, but 
has the disadvantage that it generates an overdetermined system, which increases 
the total computation time. 


4.40 Example. We illustrate this in the case considered in the Examples 4.26 and 
4.38, where we parameterized the upper half of a circle of radius one. There we 
stated that we must use the two corners of the boundary to determine the primary 
knots tf; = 0, 1, 2, and that the corresponding smoothness conditions depend on the 
approximate density @[y/t]]. The latter is expected to be discontinuous at the pri- 
mary knots. We construct a B-spline basis by adjoining secondary knots. In Example 
4.38, this was done by placing two secondary knots ¢; ; between consecutive primary 
knots, raising the total number of knots to 6. We then constructed a B-spline basis 
using cubic polynomials with continuous second-order derivatives at the secondary 
knots and no continuity at the primary knots. This basis, which consists of 12 ele- 
ments, is shown in Fig. 4.24. The primary and secondary knots are marked on the 
graph on the left in Fig. 4.26. 
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0 1 2 0 1 2 


Fig. 4.26 Left: the knot locations on 0S. Right: the collocation points. 


We can select collocation points x; € 0S to determine the 12 unknown coefficients 
Co,,i,j- Because each equation of the boundary element method has two components, 
we need at least six collocation points. We choose these points midway between 
adjacent knots. Of course, the choice is not unique, but this particular one provides 
a good position for the collocation points since it makes the integrand D|x;, y[r]] 0 
@|y[t]] ‘well behaved’ in the vicinity of x;. Our 6 collocation points are marked on 
the graph on the right in Fig. 4.26. 


4.13 Code Validation 


The four primary steps of this essential process in the construction of the code for 
our methods are 


(i) theoretical code validation; 

(ii) visualization of the code; 
(iii) validation of both small and large code parts; 
(iv) validation with known test examples. 


Each step of the code development implements a specific portion of the theory for 
solving our boundary integral equations, which means that we can often make a 
theoretical check when the step has been completed. For example, as seen above, 
computing D and P generates a large amount of symbolic code. It is known that 
these matrices must satisfy (4.6) and (4.10). Frequently, this can be verified by di- 
rect substitution and use of the symbolic simplification capabilities of Mathemat- 
ica®. It should be understood that the functions Simplify and FullSimplify 
return the primary generic result of a simplification, excluding such issues as sin- 
gular points. These functions do not locate the singularities and do not generate 
6[x,y]. Both D and P are extremely large, which makes Z, o D[x,y] and Z,o P[x, y| 
even larger. Consequently, Simplify and FullSimplify can take a substantial 
amount of time to reach the desired result. An explanation of the scope of these two 
functions was given in Remark 4.15. 

As a tule, it is desirable to have numerical validation in addition to symbolic 
theoretical validation. The best way to validate a substantial amount of numerical 
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calculations is through visualization. Sometimes, however, this is not practical, but 
in most circumstances it can be readily enacted. 


4.41 Example. Visualizing the requirement that D and P must satisfy (4.6) and 
(4.10)—that is, 


00 
Z,° Dix, y] = (( ) ? 


00 
Zxo Pix, y| = & ) , 


is straightforward. A downside of using numerical visualization is that we need to 
assign specific values to all the parameters in D and P in order to have the symbolic 
expression reduce to a numerical one, and it is always possible that such specific 
test values may be either unrepresentative or inappropriate. In our case, Fig. 4.27 
illustrates visually the condition for Z, © D[x,y] with x = {0.5, 0.5}. These graphs 
show that each of the four components is 0 except possibly near the singularity. Be- 
cause a large portion of the boundary integral equation code is numerical in nature, 
visualization becomes one of the most important aspects of the validation process. 


Fig. 4.27 Graphical representation of Z, oD[x,y], x ES. 
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It is also important to validate both small and large portions of the code. Valida- 
tion of small code segments helps ensure that the code has been created properly; 
validation of larger code segments guarantees that the individual smaller compo- 
nents are working together as required. Small code validation is often easy and 
straightforward, and can usually be accomplished by both symbolic and numeri- 
cal procedures. Every effort should be made to implement as many such validation 
tests as possible since a small segment of code, owing to its brevity, is normally well 
understood and easy to debug. Delaying debugging to larger code segments is unde- 
sirable because it is often unclear what portion of the code is causing the problem. 
However, validation of large code segments remains important, as it enables us to 
verify the overall functioning of the program. 

The validation process must finish with a number of test cases that implement the 
full and complete code. For our model, a test case should be a nontrivial solution of 
the homogeneous problem 


Zxo ulx] = (o : (4.25) 


Following a procedure similar to that used to derive the fundamental matrix of so- 
lutions, we assume that the solution we are seeking is of the form 


u(x] = adjoint [Z,] ts (4.26) 


Inserting this expression in the above equation yields 


Z, adjoint [Z,]o fee ) = Det [Z,]o ee ft) = ) 


testy [x] testy [x] 


so the two components of the solution will be constructed from the solution of the 
equation 


Det [Z,] o testy [x] = 0, 


which, as shown in Sect. 4.4, is the same as 


(A +2q) (4 53) (4 


k 
7+ 2u 3) © testa [x] = 0. 


Changing to polar form and continuing the reasoning in that section, we arrive at 


the linear combinations 
es aa, sf) =a e€SSe. {r esse. yf 
aa (04 (04 a y) 2 


where, in view of the homogeneous nature of (4.25), we have chosen the nonsingular 
function Besse1T instead of Besse1K. Once again, the code can be used to check 


. (4.27) 
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the accuracy of the above expression by verifying symbolically that 
Det [Z,] o testy [aa,ba,x] = 0. 


All the ingredients are now available to compute the desired test solution (4.26). 
Using FullSimplify and, say, the specific values 


{ai,a2} = {1,—3}, {b1,b2} = {2,4}, 


we reduce the two components ug,[1,2,—3,4,x1,x2] to 


{ (Ka +)(2x2(A +2) 


HypergeometricOF1Regularized | 1, 


k((1 ee) 


4u 
(x2a, — (1+x1)ay) +(A+2p) 


k((1+x1)?+x2? 
HypergeometricOF1 Regularized 2. (( ) | 


4u 

((1+x1 —x2)(1+x1+x2)a, +2(1+x1)x2a)) —2(1+x1)y 
eae) 

4(A + 2) 


k 
HypergeometricOF 1 Regularized f ; (( 


((1+x1)b) +x2b2) + 


k((1+x1)? +x2? 
HypergeometricOF1 Regularized [2 (( ) | 


4(A +2) 
((1+x1 —x2)(1+x1+x2)b, +2(1+x1)x2bp))) / 
(2((1+x1)? +x2*)u(A+2y)), 


(K(A + w)(2(1+x1)(A + 2p) 


k((1+x1)? +x2? 
HypergeometrcOF Regularize| (( ie 


4u 
(—x2a) + (1+x1)ar) +(A +2) 


k((1+x1)? +x2? 
HypergeometricOF | Regularized [2 ( ) | 


4u 
(2(1+x1)x2a; + (—(1+x1)? +x2”)az) — 2x2 
k((1 ee) 
4(A + 2u) 


HypergeometricOF 1 Regularized | 1, 
((1-+x1)b) +x2b2) + 


k((1+x1)? +x2’) 
HypergeometricOF1Regularized | 2, 


4(A + 2u) 
(2(1 +x1)x2b; + (— (1+x1)? +x27)b2)))/ 


(2((1+x1)? +x27)u(a +2u)). 
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This test solution can be verified symbolically by applying FullSimplify 
to check that it satisfies (4.25). The verification can also be done graphically, as 
illustrated in Fig. 4.28. 


1.x 107! 


-1.x107! 


Fig. 4.28 Graphs of the evaluation of the components of Z, oul[a1, Bi, 02, Bo, x]. 


Full verification of the code can be achieved by considering a wide range of 
particular test cases. 


4.42 Example. Consider the homogeneous equation in S with a Dirichlet condition 
on OS; that is, 


The data function Y is constructed from a known test solution u[aj,b1,a2,b2,x] 
restricted to the boundary. We use the test solution generated with the parameter 
choice 


A>1, wo2, k-3, 


aml, Bo2, m>-3, b-4. 


Our test domain S is the half disc featured in Example 4.26 (see Figure 4.18), with 
the same parameterization (4.21) of the boundary 0S. The graphs of the two com- 
ponents of the test solution u in S are shown in Fig. 4.29. 

The two components of Y on dS can be constructed directly from u[x] by re- 
stricting x to the boundary: 


Px) =ulx], xe OS. 


The graphs of the components of this function are displayed in Fig. 4.30. 
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—100 


Fig. 4.30 The components of 7. 


We solve our Dirichlet problem by means of the direct method described in Sect. 
3.2, which yields the boundary equation (2s ) (see Example 4.24). This equation 
is coded in the form (4.19), where the unknown density @ represents Tu. Once @ is 
determined, we find the solution in § from the representation formula (4.20). 

We can use our test solution wu to construct the exact function Tu by evaluating 


(Tu) ik] = Fouls] = To ee = eae , xedS. 


up [x] Tuo |x] 
The graphs of both components of Tu are shown in Fig. 4.31. 


2 


Fig. 4.31 The components of Tu. 
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An approximation @ of @ = Tu can be calculated by implementing a boundary 
element method of our choosing. As we mentioned earlier, our preferred choice is 
the collocation method in conjunction with B-spline basis elements b; ;. The func- 
tion @ will have discontinuities at the two domain corners, and the approximating 
B-spline basis will of necessity reflect this feature. For our problem, we use a piece- 
wise cubic spline that is twice continuously differentiable at the secondary knots 
and discontinuous at the corner knot locations t = 0, 1, 2. In terms of the selected 
parameterization, the primary and secondary knot locations and the corresponding 
smoothness conditions are specified by the set 

{0, 0, 0,0, 4, 5,3) 1,1, 1,1, 2, 3, 7,2, 2, 2, 2}. 

These locations, marked on the graph on the left in Fig. 4.32, determine 14 functions 
b;,; for each component of @. The graphs of these functions are shown in parametric 
form in Fig. 4.33. 


0 1 2 0 1 2 


Fig. 4.32 Left: the knot locations on 0S. Right: the collocation points. 


We need to choose a suitable number of collocation points to be able to compute 
a solution @. This can be done in several different ways. But we must remember 
that a collocation point cannot be placed at any of the corners t = 0, 1, 2, where @ 
is expected to be discontinuous. The 14 parametric values for t used as collocation 
points and displayed on the graph on the right in Fig. 4.32 are 


{2 3.5 1 9 I 13 15 17 19 3 23 25 27 
14> 14? 14> 2) 14> 14> 14) 14> 14> 14? D> Ta 14> 14S 


The same 14 collocation points are used for both components of @ to produce 28 
equations, from which we compute the coefficients cg; ; in the representation 


Olt] = {Grlx{e]], Golxle]]} 


= {LDevamislhTDeassP tl} 0<t<tMax. 
J J 


i i 


The graphs of the two components of @ are shown in Fig. 4.34. This approximate 
solution agrees very well with the exact solution @ = Tu (see Fig. 4.31). 
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Fig. 4.33 The 14 B-spline basis functions b,, ;[t], 0<t < 2. 
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Fig. 4.34 The components of @. 
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To get a better understanding of the discrepancy between @ and @Q, in Fig.4.35 
we have graphed the difference @ — @ relative to the maximum absolute value of @ 
on OS. 


Fig. 4.35 The components of the relative error in @ (Cartesian coordinates). 


Since Fig. 4.35 is somewhat difficult to interpret, in Fig. 4.36 we have graphed 
the same function in parametric form for 0 < t < 2. These new graphs indicate that 
the relative error in the two components of @ is only about 0.2%, except possibly 
near the corners of 0S. 
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Fig. 4.36 The components of the relative error in @ (parametric form). 
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We can now use @ and # together with the representation formula 
ie] = $ Dey] GH] at, f Phe,yloubl dh, xe, 
r r 


to compute an approximation # of the exact solution u in S. The graphs of the two 
components of i are displayed in Fig. 4.37. Comparing these graphs with those of 
the components of u in Fig. 4.29, we conclude that the difference between # and u 
cannot be discerned visually. 


—100 


N 
2 


Fig. 4.37 The components of a in S. 


4.43 Remark. Example 4.42 and many others presented in the next chapter for all 
types of boundary conditions and boundary integral methods validate the computer 
code of our Mathematica® program. 


Chapter 5 
Computational Examples 


5.1 Preliminaries 


5.1.1 Test Solution 


As established in the earlier part of the book, the equilibrium system governing our 
mathematical model is 


Zyoulx] =0, xES, (5.1) 


where u is the displacement vector function and Z,. is the operator defined by (4.1). 
The examples discussed in this chapter illustrate the numerical implementation of 
the direct and classical indirect methods for various boundary value problems. Since 
it is important to know how accurate our results are, in a majority of cases we make 
use, for comparison purposes, of a test solution of system (5.1). Specifically, using 
(4.26) and (4.27) with parameters 


As1, p32, k33, m=1, b=2, m=-3, b=4 (5.2) 


in conjunction with the Mathematica® function FullSimplify, we construct 
the particular solution 


u|x1, x2] 


= { i ((2 x2*HypergeometricOF1 (1, 3 ((1 + x1)? +x2’)| 


+ (1+xl—x2)(1+xl14 x2)HypergeometricOF I [2, j ((1+x1)?+x2*)]) / 


((1+x1)* +x2?) 
+4(- (2c1 +x1)PBesseut|2,y/$/(1 4x1)? 422"])) / (14x)? +902) 


— HypergeometricOF 1 Regularized (2, a (+x 4x2? ) i) 
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= 8 (1 +x1)x2 HypergeometricOF1Regularized [3,4 (+ x1)? +x2? )] 


+7 (1 +x1)x2 HypergeometricOF1Regularized (3, 3 ((1+x 1)? + x2°)| ) ; 


3 (s ( — 2 x2°>HypergeometricOF1Regularized [1,4 (+ xl) + x2”)| 
+ (—(1+x1)*+x27) 
HypergeometricOF 1Regularized{2, 3 30 (( +x1) + x2’) | ) / 
(5((1+x1)?+x24) ) 
=o 6 (1 +x1)x2 HypergeometricOF1Regularized [3,4 my (A+ xl)? + x2°)| 
- 3 (1+x1)x2 HypergeometricOF 1Regularized|[3, 2 (( (1+x1)?+ x2°)| 
- 3 (4 HypergeometricOF 1Regularized (2, ; 2((1 +x1)?+ reall 
+3(1+x1)*HypergeometricOF1Regularized [3, 3 (+x 4 x2?)))) \ : 
(5.3) 


After simplification, the corresponding expression for (Tu)[x], constructed with 
(4.2) and the same parameter values (5.2), is 


(Tu) [x1, x2] 


= {9(so((1 + x1)? +6(1 +x1)?x2 —3(1+x1)x2? — 2x2?) 


HypergeometricOF1 | 1, A ((1+x1) +x2’)] 
+ 100((1-+x1)? —9(1 +x1)?x2— 3(1 +x1)x2? +3 x23) 
HypergeometricOF1 |1, 3 ((1 +x1)? +x2°)| 
2(5(1 + x1)3(11 +3 x1(2-+x1)) +30(1-+x1)2(9-+x1(2-+x1))x2 
+ 6(1-+x1)(—17+3x1(2+x1))x2? +4(—114+9x1(2+x1))x2? 
+3(1+x1)x2*+6 x25) HypergeometricOF [2, 34 ((1+x1)* +x2’) | 


25( 4(1+x1)? +9(1 +x1)?(5+x1(2+x1))x2+43(1+x1) 


(5-+x1(2-+x1))x2? +3(—143x1(2+x1))x23 +3(14 x1)x2*) 
HypergeometricOF 2, 3 ((1+x1)* +x2?)]) vi [xt,x2] / 


(50((1+x1)? +x2”)”) 


9( 160(2(1 +x1)? —3(1 +x1)2x2— 6(1 + x1)x2? + x23) 
HypergeometricOF 1 Regularized|1, a (A+ x1)? + x2°)| 
+ 200(3(1 -+x1)? +3(1+x1)?x2—9(1 +x1)x2? — x2?) 
3 ((1+x1)? +x2’)| 
16( 20(1 +x1)3 +3(1-+x1)2(11-+x1(2+x1))x2+6(1+x1) 


HypergeometricOF1 Regularized [1, 
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(11 +x1(2+x1))x2? + (—743 x1(2-+x1))x23 +6(1 + x1)x2*) 


HypergeometricOF1Regularized|2, x (+x PP +x2? ) 


25(3(1 + x1)3(11 +3 x1(2-+x1)) +3(1 +x1)?(9+x1(2-+x1))x2 


72(1 + x1)x22 — 8 x23 —9(1+x1)x24—3 x2°) 
HypergeometricOF 1Regularized|[2, 2 (( +x1)?+ x2°)| ) v2fx 1x2] / 


(100((1 + x1)? +x22)”), 


9( 160(2(1 +x1)3 — 3(1 +x1)2x2 — 6(1 + x1)x2? + x23) 
HypergeometricOF1Regularized|1, * i; +x1)?+x2? ) 
+ 200(3(1-+x1)? +3(1 +x1)?x2 — 9(1 +x1)x2? — x23) 
HypergeometricOF1 Regularized (1, 3 (1+ xl)? + x2’)| 
16( 20(1-+x1)? +3(1 +x1)?(11 +x1(2+x1))x2+6(1+x1) 


(11+ x1(2-+x1))x2? + (—7+3 x1(2-+x1))x2* +6(1 +x1)x2*) 


HypergeometricOF1Regularized (2, * ((1+x1)? +x2’)| 
25(3(1 +x1)3 (11-43 x1(2+x1)) +3(1 +x1)2(9-+x1(2 +x1))x2 


72(1 +x1)x2? — 8 x23 —9(1-+x1)x24 3x28) 
HypergeometricOF | Regularized (2, 3 ((1 +x1)? +x27)| ) v1[x1,x2] / 


(100((1 +x1)? +x22)”) 


9(80((1 + x1)? +6(1 +x1)2x2—3(1-+x1)x2? —2 x23) 
HypergeometricOF 1 [1, 3 Le +x2’)] 
+ 100((1-+x1)? —9(1 +x1)?x2 — 3(1 +x1)x2? +3 x23) 
1+x1)?+x2’)] 
2+x1)) +6(1+x1)?(—39+x1(2+x1))x2 


HypergeometricOF1 [1,3 2 ( 
+2((I +x1)3(—37+3x1 


( 
( 

+ 6(1+x1)(23 +3 x1(24 et 
)) 


+ 4(2949 x1(2-+x1))x23 + 15(1 + x1)x2* +30 x2°) 
HypergeometricOF1 [2, # (( (( 1+x1) ee 


25( 4(1+x1)? +9(1 +x1)?(5+x1(2+x1))x2+3(1+x1) 


(S+x1(2-+x1))x2? +3(—1+3x1(2-+x1))x2* +3(1 +x1)x24) 


HypergeometricOF I [2, 3 ((1+x1)* +x2?)]) v2[x1,x2] / 
(so((1 +x1?-+x22)?) b (5.4) 


where v1[x1,x2] and v2[x1,x2] are the components of the unit normal to 0S. 
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5.1 Remark. The expressions of both u[x] and (Tu) [x] are constructed with the func- 
tions HypergeometricOF1 and Hypergeometric0OF1lRegularized. The 
former is the confluent hypergeometric function oF\(a;z) and the latter is the reg- 
ularized confluent hypergeometric function 9F{(a;z)/T (a) (see the Mathematica® 
documentation). 


5.1.2 Computational Accuracy 


In Chap. 4, we explained in great detail that the evaluation of singular integrals 
requires a significant amount of computational attention. As a consequence, the en- 
tire code for the examples presented here is written using exact arithmetic. Floating 
point machine accuracy is usually limited to about 16 digits. However, when han- 
dling the singularities occurring in our problems, we often need to compute with 
25 or more significant digits. This places a severe restriction on how the code is 
written. In Mathematica®, the accuracy of an expression is always dictated by the 
accuracy of its least accurate term. Consequently, it is impossible for the entire code 
to contain floating point numbers with a predefined fixed accuracy. 

There is a difference between representational accuracy and computational accu- 
racy. Suppose that the number 7 is to be entered with a representational accuracy 
of five digits. The floating point number 3.1415 can be entered with full represen- 
tational accuracy as 31415/10000. On the other hand, suppose that our computation 
of a difficult expression would result in a loss of 10 digits. Then, in order to pre- 
serve the five digits of representational accuracy for our number, the computational 
accuracy must be increased to 15 digits. 

The examples in this chapter have been computed with a representational accu- 
racy of 10 digits. For the “better behaved’ singularities encountered in them, there 
is normally a loss of less than five computational digits; hence, for these examples, 
16-digit machine accuracy is sufficient. However, in some ‘ill behaved’ cases, a loss 
of 10 digits or more is possible, so, since here machine accuracy is not adequate, the 
computational accuracy must be increased to at least 20 digits. 

In Mathematica®, the representational accuracy of a computation for a quan- 
tity x is controlled by PrecisionGoal and AccuracyGoal. The distinction 
between the two is determined by how the error is bounded; that is, 


|error| ss 1Q~AccuracyGoal te Eel Ftaensort 


Computational accuracy is specified by the setting for WorkingPrecision, 
which is defined as the number of digits used in the computation. Normally, the 
WorkingPrecision must exceed the values of both the PrecisionGoal and 
AccuracyGoal by 5 to 10 digits; in certain circumstances, this number may be 
even larger. 

The above discussion refers to controlling the loss of accuracy for a single com- 
putational step. But the final answer is usually the result of a large number of such 
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steps, each involving a loss of accuracy. This loss will accumulate to produce an 
answer with less than the desired representational accuracy. So, for example, if five 
digits of accuracy are desired in the final answer of a procedure where six dig- 
its are lost due to the cumulative computation and seven digits may be lost for a 
single difficult computation, then the WorkingPrecision should be set to at 
least 18 digits when performing that difficult computation. We emphasize that the 
largest part of the program produces only a small loss of accuracy, which means that 
only a limited portion of the code requires a high value of WorkingPrecision. 
The WorkingPrecision, PrecisionGoal, and AccuracyGoal have lim- 
itations, and their proper use requires careful consideration, numerical validation, 
and continual vigilance. 


5.2 Remark. The examples discussed in the rest of this chapter share a common, 
streamlined template. To avoid forcing the reader to flip back through the pages 
in search of various quoted formulas, we have made each example self-contained 
by concentrating the necessary (analytic, numerical, and graphical) information in 
each individual case and adopting a generic economical style of presentation. After 
considering various other options, we felt that a certain amount of repetitiveness was 
a small price to pay to have the examples easily and independently accessible. 


5.2 Dirichlet Problem in an Ellipse: Piecewise Linear Spline 
5.2.1 Summary 


This example illustrates the application of the direct method to the Dirichlet prob- 
lem in a convex symmetric domain with a smooth boundary curve. A known test 
solution is used to calculate the corresponding boundary data function. The collo- 
cation method is then implemented, with piecewise linear splines on equally spaced 
collocation points, to construct an approximation of the original solution from its 
boundary values. The results are validated through a comparison between the exact 
and approximate solutions. 


5.2.2 Problem Statement 
Domain boundary. 0S is the ellipse parameterized by 
x1 [tf] =1+Cos[xt], x2[r] = 5Sin[ar]+5, O<1<2. (5.5) 


Its graph can be seen in Fig. 5.6. 


Governing equations. The displacement u is the solution of the boundary value 
problem 


Zyoulx] =0, xES, 
ulx}= Pix], xe. 
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Test solution. The two components of the test solution u given by (5.3) are graphed 
in Fig. 5.1. 


Fig. 5.1 The components of u[x], x € S. 


Boundary data function. The function Y is computed from the test solution wu as 


Pix] =ulx], xeS. 


The graphs of its two components are displayed in Figs. 5.2 and 5.3, in Cartesian 
coordinates and parametric form, respectively. 


Fig. 5.2 The components of Y [x] (Cartesian coordinates). 
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Fig. 5.3 The components of Y[x{¢]] (parametric form). 
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5.2.3 Solution Procedure 


Method. We apply the direct method, which reduces the problem to the boundary 
integral equation 


Vo(p) = (Wo+51)A onads, 


coded as 


$ Dix, y|o oly] dT, = f Plx,ylo Plyldl +3 Pil. 


q Weak(x] 


(5.6) 
TCPVv is] 


This equation is solved numerically to obtain an approximate density @. In turn, @ 
is used in the representation 


= $ Dix sogbldty— $ Plx.xjoPhldh, xes 
FE Tr 


to generate an approximate solution @ in S. 


(5.7) 


Density function. The exact density @ is the boundary stress vector Tu computed 


from the test solution u and given by (5.4); its components are graphed in Figs. 5.4 
and 5.5. 
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Fig. 5.5 The components of @[x{t]] (parametric form). 
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Numerical approximation. We compute @ by the collocation method with a 
B-spline basis consisting of elements b;,;. Then the approximate density is sought 
in the form 


Lhe, jij] 
ij 


5.8 
LDXc24,jb1,,lt] J’ ee) 
bey 


$lalt]] = 


where the numerical coefficients cy ;,; are determined by substituting (5.8) in equal- 
ity (5.6). 

Since @ is continuous on the boundary, we use a piecewise linear spline with 
continuity at all knot locations. In terms of the parameterization (5.5), the primary 
and secondary knots and the smoothness at their locations are specified as the set 


{0 123 4 5 6 7 8 9 10 1 
10? 12% 11% 11% 11% 112 2h? 11s Wis 179 


12 13 14 15 16 17 18 19 20 21 2} 
127 117 117 12° 12° 12? 11? 11? 11? 11? . 


The knots are marked on the graph on the left in Fig. 5.6. 
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Fig. 5.6 Left: the knot locations on 0S. Right: the collocation points. 


The 22 knots determine 22 functions b;,; for each component of @, generating a 
total of 44 basis functions. The graphs of the b;,; in parametric form are shown in 
Fig. 5.7. 

To compute ©, we must choose an appropriate number of collocation points, 
which can be placed anywhere because D, P, @, and x are smooth on 0S. The set of 
values of the parameter ¢ at our selection is 


{0 123 4 5 6 7 8 9 10 
2127 11? 127 11? 12? 11? 11? 11? 11? 11? 


1. 12. 13 14 15 16 17 18 19 20 21) 
> TT? Tr? TW? Wr? TW? TW? TW? TW? WW? TI 
These values are marked on the graph on the right in Fig. 5.6. 
The same 22 points are used for both components of @, yielding a system of 44 
constraining equations that enables us to compute the coefficients cq ; ; in (5.8). The 
44 x 44 matrix of this system has a condition number of 24. 
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Fig. 5.7 The 22 B-spline basis functions 0j,;[t], 0<t < 2. 
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5.2.4 Solution 


Approximate density. The two components of @ are graphed in Figs. 5.8 and 5.9. 
The latter also shows the spline knot locations. It can be seen that @ agrees very 


well with the exact density @ in Fig. 5.4. 


100; 
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Fig. 5.9 The components of [x[t]] (parametric form) and the knot locations. 


Approximate solution. We use @ and # in (5.7) to compute an approximation a 
to the exact solution u in S. Its two components are graphed in Fig. 5.10. 


—100 
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Fig. 5.10 The components of ai[x], x € S. 
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Error analysis. To visualize the difference between the approximate and exact den- 
sities, in Fig. 5.11 we graphed @ — @ relative to the maximum absolute value of @ 
on oS. 

The same was done in parametric form in Fig. 5.12, which indicates that the 
relative error in the two components of @ is only about 1%, except possibly at the 
extreme right of the boundary, where it can exceed 3%. 


1/2 1 3/2 
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Fig. 5.12 The components of the relative error in @[x[t]] (parametric form). 


Figure 5.11 displays a larger error at the extreme left and right segments of the 
boundary, where 0S has the largest curvature and where, therefore, the integrand is 
changing most rapidly. Also, Fig. 5.11 is constructed from the error relative to the 
maximum absolute value of @ on OS, which creates distortion and leads to a larger 
relative error at the right side of the boundary, where @ attains its maximum absolute 
value. 


5.3 Remarks. (i) The B-spline basis functions had to be modified to produce a basis 
with periodic continuity at the knot locations t = 0, 2. The first function in Fig. 5.7 
illustrates the required modification. 
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(11) We notice that the equally spaced parametric knot locations and collocation 
points are, to all intents and purposes, also equally spaced geometrically over the 
boundary (see Fig. 5.6). Consequently, the error is expected to be essentially the 
same along the entire boundary contour (see Fig. 5.11). In general, the position of 
the knots and collocation points should be determined by the boundary 0S rather 
than by its parameterization. 


(iii) The set of coefficients {cy ij, ¢2,,;} in (5.8), computed with the B-spline 
basis shown in Fig. 5.7, is 


{{40.1, —44.5, 83.8, 87.3, —75.4, 59.4, —44.8, —33.4, 
25.1, -18.9, —12.9, —4.39, 5.41, 13.4, 21.2, 30.8, 43.8, 61.1, 
82.9, 106., 121., 109.}, {-173., —194., —170., —136., 
0.953, 14.5, 21.3, 21.7, 20.6, 19.4, 17.8, —104., —77.4, —55.8, 
— 38.8, —25., —12.4, 14.5, 7.08, —9.24, —42.3, —102.}} 


5.3 Dirichlet Problem in an Ellipse: Piecewise Constant Spline 


5.3.1 Summary 


The difference between this example and the preceding one lies in the nature of 
the functions that form the B-spline basis. Using a piecewise constant spline, we 
construct a family of basis functions with discontinuities at the knot locations and 
comment on how to make an appropriate choice of (equally spaced) collocation 
points. 


5.3.2 Problem Statement 


Domain boundary. QS is the ellipse parameterized by 


x1 [t] =1+Cos[zt], 


; 0<1< 2, (5.9) 
x2[¢] = 5 Sin[xt] +5, 


Its graph can be seen in Fig. 5.18. 


Governing equations. The displacement wu is the solution of the boundary value 
problem 


Zyoulx] =0, xES, 
ulx]}= Pix], xe. 
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Test solution. The two components of the test solution u given by (5.3) are graphed 
in Fig. 5.13. 
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Fig. 5.13 The components of u[x], x € S. 


Boundary data function. The function Y is computed as 
Px] =ulx], xe Os. 


The graphs of its two components are displayed in Figs. 5.14 and 5.15, in Cartesian 
coordinates and parametric form, respectively. 


Fig. 5.14 The components of [x] (Cartesian coordinates). 
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Fig. 5.15 The components of [x(t]] (parametric form). 


106 


5 Computational Examples 


5.3.3 Solution Procedure 


Method. We apply the direct method, which reduces the problem to the boundary 
integral equation 


Vo(~) =(Wo+51)A ond, 
coded as 


f Pkerlooblah= f PhrloPhldh +3 Abi. 
] 


TWeak(x 


(5.10) 
TCPVv ix] 


This equation is solved numerically to obtain an approximate density @. In turn, @ 
is used in the representation 


ul] = $ Dislohldty— $ PixyloPbldh, xeS, 
Tr Tr 


to generate an approximate solution @ in S. 


(5.11) 


Density function. The exact density @ is the boundary stress vector Tu computed 


from the test solution u and given by (5.4); its components are graphed in Figs. 5.16 
and 5.17. 
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Fig. 5.16 The components of @[x] (Cartesian coordinates). 
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Fig. 5.17 The components of ¢[x[t]] (parametric form). 
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Numerical approximation. We compute @ by the collocation method with a 
B-spline basis of elements b;,;. Then the approximate density is sought in the form 


Ldc1 i, /bi,;{t] 


[x(t] = Se (5.12) 
ij 


where the numerical coefficients cq ;,; are determined by substituting (5.12) in 
(5.10). 

We use a piecewise constant spline for which, in terms of the parameterization 
(5.9), the knot locations and the corresponding smoothness conditions are specified 
as the set 


123 4 5 6 8 9 10 11 12 13 
{0, 5, 3,3, 9,3,9,1, 5,2, 2,4, 2, 3,2}. 


The knots are marked on the graph on the left in Fig. 5.18. 
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Fig. 5.18 Left: the knot locations on 0S. Right: the collocation points. 


The 14 knots determine 14 functions b;,; for each component of @, generating a 
total of 28 basis functions. The graphs of the b;,; in parametric form are shown in 
Fig. 5.19. 

To compute @, we must choose an appropriate number of collocation points. 
Since D, P, @, and x are smooth on 0S, such points can be placed anywhere except 
at the knots, where the b; ; are discontinuous. Our choice is to position them midway 
between the knots, which means that the 14 values of t selected as collocation points 
are the elements of the set 


{2 3.5 1 9 I 13 15 17 19 3 23 25 O74: 
14? 14> 14? 2? 14? 14> 14? 14? 14> 14? 2? 14> 14> 14 
These values are marked on the graph on the right in Fig. 5.18. 

The same 14 points are used for both components of the approximate density @, 
yielding a system of 28 constraining equations that enables us to compute the nu- 
merical coefficients cg j,; in (5.12). The 28 x 28 matrix of this system has a condition 
number of 11. 
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Fig. 5.19 The 14 B-spline basis functions bj, ;[t], 0<t <2. 


5.3.4 Solution 


Approximate density. The components of @ are graphed in Figs. 5.20 and 5.21. 


Refined numerical approximation. From Fig. 5.20 we see that the relative error 
in the components of @ exceeds 20%, which means that we need a finer spline knot 
spacing. To achieve this, we consider the new knot locations 


{0,2 23 4167 8 9 211 12 13 14 3 16 
> 35> 25> 25> 259 5% 25> 25? 25> 25> 59 259 259 259 257 5? 25? 


17 18 19 4 21 22 23 24 7 26 27 28 29 6 31 32 33 
257 257 2579 5% 25? 25? 25? 25? ~? 257% 25? 25? 25? 5? 25? 25? 25° 


34 7 36 37 38 39 8 41 Fey Bos a5) 2) Geo aes ger ge 2h 
25332 25% 257 25:9 25% 5225.7 25:2 25:2 25% 5.2 259257 257 25? . 
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Fig. 5.21 The components of @[x[t]] (parametric form). 


These knots are marked on the graph on the left in Fig. 5.22. 


0 I 2 0 1 2 
Fig. 5.22 Left: the knot locations on 0S. Right: the collocation points. 
Our chosen 50 collocation points necessary to compute @ are marked on the 


graph on the right in Fig. 5.22. The corresponding values of the parameter ¢ that 
generate them are the elements of the set 
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{ 5s Sys iO SD? Sh 0) 4B? Ta? a> Be) aa? SO D> Sar ab» 2o> Be 
50’ 50’ 10’ 50° 50’ 50° 50° 10’ 50’ 50’ 50? 50? 2? 50° 50’ 50’ 50° 


7 37 39 41 43 9 47 49 51 53 Il 57 59 61 63 13 67 
10’ 50’ 50’ 50° 50’ 10’ 50’ 50’ 50’ 50’ 10’ 50’ 50’ 50’ 50’ 10° 50° 


69 71 73 3 77 79 81 83 17 87 89 91 93 19 97 99 
50°? 50? 50° 2? 50° 50° 50? 50? 10°? 50? 50? 50’ 50” 10° 50° 50 J° 


Using the same 50 collocation points for both components of @, we arrive at 
a system of 100 constraining equations for the coefficients cg; ; in (5.12). The 
100 x 100 matrix of this system has a condition number of 44. 


Refined approximate density. Figure 5.23 shows the graphs of the two compo- 
nents of @ computed with the new setup. 
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Fig. 5.23 The components of @[x] (Cartesian coordinates). 


The components of @ and @ are displayed together in Fig. 5.24. These graphs 
show good agreement between the two functions at the collocation points, which 
were placed midway between the new spline knots. 
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Fig. 5.24 The components of @[x[t]] and @[x[t]] (parametric form). 
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Approximate solution. We use the refined version of @ and Y in (5.11) to compute 
an approximation @ to the exact solution u in S. The graphs of the two components 
of # are shown in Fig. 5.25. 


-100 


Fig. 5.25 The components of [x], x € S. 


Error analysis. To visualize the difference between the approximate and exact 
densities, in Fig. 5.26 we graphed @ — @ relative to the maximum absolute value of 
@ on OS. The same was done in parametric form in Fig. 5.27, which indicates that 
the relative error in the two components of @ is only about 5%, except possibly at 
the extreme left and right sides of the boundary, where it can exceed that value. 


Fig. 5.26 The components of the relative error in @[x] (Cartesian coordinates). 
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Fig. 5.27 The components of the relative error in @[x[¢]] (parametric form). 
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5.4 Remarks. (i) We notice that the equally spaced parametric knot locations and 
collocation points are, to all intents and purposes, also equally spaced geometrically 
over the boundary (see Fig. 5.18). Consequently, the error is expected to be essen- 
tially the same along the entire boundary contour (see Fig. 5.23). In general, the 
position of the knots and collocation points should be determined by the boundary 
OS rather than by its parameterization. 

(ii) The set of coefficients {c1 jj, c2;,;} in (5.12), computed with the B-spline 
basis shown in Fig. 5.19, is 


(i975; 877,766, -519,—32.8,-21, 
— 11.2, 3.31, 16.1, 30.1, 51.9, 84.7, 118., 97.1}, 
{—194., —158., —106., —66., —37.8, —16.9, 4.19, 
20.5, 21.5, 19.6, 16.5, 6.26, —28.3, -119.}}. 


The corresponding set for the B-spline basis constructed with the refined selection 
of knots and collocation points displayed in Fig. 5.22 is 


{{19.8, —17.2, —47.5, —68.5, —80.9, —86.4, —86.9, 
83.9, —78.8, —72.4, —65.5, 58.6, —52., —45.8, —40.3, 
45.4, —51.2, —27.45, —243, 214, —18.8, —16.2, 
13.4, —10.2, —6.4, —2.14, 2.15, 6.18, 9.85, 13.2, 
16.6, 20.2, 24., 28.3, 33.1, 38.7, 45., 52.2, 60.3, 69.2, 
78.9, 89., 99.1, 108., 115., 119., 118., 108., 88.5, 57.5}, 
{—180., —191., —190., —182., —170., —156., —141., 
126... =112,,—99.4, =87 4, —76.4, 66.3, —57.4, 
49.2418, —35:2. —29., —23.3,.-17.7, 12.4, 
—644,—0,45,5.7, 11:5, 16.3; 19.5, 21.2, 91,7, 21.6, 
21.3, 20.8, 20.3, 19.7, 19.1, 18.4, 17.5, 16.3, 14.5, 12., 8.36, 
3.2, —3.99, —13.9, —27.3, —45.1, —68.1, —96.3, —127., —158.}}. 


5.4 Dirichlet Problem in an Asymmetric Domain: Piecewise 
Linear Spline 


5.4.1 Summary 


This case is similar to the one discussed in Sect. 5.2, with the difference that here 
the direct method is applied to the Dirichlet problem in an asymmetric domain with 
a smooth boundary. 
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113 
5.4.2 Problem Statement 
Domain boundary. QS is the curve parameterized by 
We=1(7 ACos|[2t 
x1[f] 0 (7 + 8Cos|zt] + 4Cos| ‘i }), 0<1<2. (5.13) 
x2[t] = <p (25 + 2Cos[3z1] + 20Sin[zz]), 


Its graph can be seen in Fig. 5.33. 


Governing equations. The displacement u is the solution of the boundary value 
problem 


Zyoulx] =0, xES, 
ulx] = Pix], xeas. 


Test solution. The two components of the test solution u given by (5.3) are graphed 
in Fig. 5.28. 
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Fig. 5.28 The components of u[x], x € S. 


Boundary data function. The function Y is computed from the test solution u as 
Px} =ulx], xe OS. 


The graphs of its two components are displayed in Figs. 5.29 and 5.30, in Cartesian 
coordinates and parametric form, respectively. 
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Fig. 5.29 The components of Y [x] (Cartesian coordinates). 
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—100 
Fig. 5.30 The components of Y[x[t]] (parametric form) 


5.4.3 Solution Procedure 


Method. We apply the direct method, which reduces the problem to the boundary 
integral equation 


Vo(~) =(Wo+51)F_ onds, 
coded as 


¢ Dkx.sJoeblar; = 
] 


TWeak(x 


$ Plx,ylo Aly|diy +4 Ply. 


(5.14) 
Topv ix] 


This equation is solved numerically to obtain an approximate density @. In turn, @ 
is used in the representation 


ul] = $ Di sloohldty— $ PixyloPbldh, xes 
Tr Tr 


(5.15) 


to generate an approximate solution [x], x € S. 


Density function. The exact density @ is the boundary stress vector Tu computed 
from u and given by (5.4); its components are graphed in Figs. 5.31 and 5.32. 


Fig. 5.31 The components of @[x] (Cartesian coordinates). 
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Fig. 5.32 The components of @[x[t]] (parametric form). 


Numerical approximation. We compute @ by the collocation method with a 
B-spline basis of elements 5; ;. Then the approximate density is sought in the form 


; LL C14 iii [] 
@|x{t]] = Dea ibij (¢] ’ (5.16) 


ij 


where the numerical coefficients cq ;,; are determined by substituting (5.16) in 
(5.14). 

Since @ is continuous on the boundary, we choose a piecewise linear spline with 
continuity at all knot locations. In terms of the parameterization (5.13), the primary 
and secondary knots and the smoothness at their locations are specified as the set 


fo, 4,4,3,2,5,3,4,4,4,5,4,5 81 
9 14> 7s 14) 79 14> 79 2 7s 14> 7 14> 7) 149 1 
15 8 17 9 19 10 3 I 23 12 25 1327, 9} 
149.7) 14> 7) 140 79 29 J 140 7? 149 7 14> 


1 1 
~o—» Om 
* ~——-0-—__» ' eo SO . 
. a ~e, Q gees 
. > S , ©) 
© “a ( 
a = aa — 
9.» -@— OW 0 Or 
0 0 
0 1 2 0 1 2 


Fig. 5.33 Left: the knot locations on 0S. Right: the collocation points. 


The 28 knots determine 28 functions b;,; for each component of @, generating a 
total of 56 basis functions. The graphs of the 5; ; in parametric form are shown in 
Fig. 5.34. 
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Fig. 5.34 The 28 B-spline basis functions b,, ;[t], 0<t < 2. 
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To compute ©, we must choose an appropriate number of collocation points, 
which can be placed anywhere because D, P, @, and x are smooth on 0S. The set of 
values of the parameter ¢ at our selection is 


4? 7? 14? 7? 14? 7? 22 7? 14? 7? 142 7° 14? 
8 17 9 19 10 3 Ll 23 12 25 13 27 
791497? 149-792? 7? 14? 7? 14° 7? 4S 


1 1 3 2 5 3 14 9 5 11 6 13 
{0, i4 I, 


15 
14? 


These values are marked on the graph on the right in Fig. 5.33. 

The same 28 points are used for both components of @[x[1]], yielding a system of 
56 constraining equations that enables us to compute the coefficients cq ;,; in (5.16). 
The 56 x 56 matrix of this system has a condition number of 74. 


5.4.4 Solution 


Approximate density. The two components of @ are graphed in Figs. 5.35 and 
5.36. The latter also shows the spline knot locations. It can be seen that @ agrees 
very well with the exact density @ in Fig. 5.43. 
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Fig. 5.36 The components of @[x[t]] (parametric form) and the knot locations. 


118 5 Computational Examples 


Approximate solution. We use @ and # in (5.15) to compute an approximation 
ii to the exact solution u in S. The graphs of the two components of @# are shown in 
Fig. 5.37. 
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Fig. 5.37 The components of ai[x], x € S. 


Error analysis. To visualize the difference between the approximate and exact 
densities, in Fig. 5.38 we graphed @ — @ relative to the maximum absolute value of 
@ on OS. The same was done in parametric form in Fig. 5.39, which indicates that 
the relative error in the two components of @ is only about 1.5%, except possibly 
at the extreme left and right sides of the boundary, where it can exceed 5%. This 
anomaly is explained by the fact that 0S has the largest curvature on those segments, 
inducing the fastest variation in the integrand and, therefore, a significant increase 
in the error. Additionally, Figs. 5.38 and 5.39 are constructed from the error relative 
to the maximum absolute value of g] on 0S, which creates distortion and leads to a 
larger relative error on the extreme right segment of the boundary, where @ attains 
its maximum absolute value. 
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Fig. 5.38 The components of the relative error in @[x] (Cartesian coordinates). 
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Fig. 5.39 The components of the relative error in @[x{¢]] (parametric form). 


5.5 Remarks. (i) The B-spline basis functions had to be modified to produce a basis 
with periodic continuity at the knot locations t = 0, 2. The first function in Fig. 5.34 
illustrates the required modification. 


(11) We notice that the equally spaced parametric knot locations and collocation 
points are not equally spaced geometrically over 0S (see Figs. 5.33). Hence, the 
error is expected to be much larger on the right side of the boundary, where these 
points are farther apart. 


(iii) The set of coefficients {c1,j,;, c2,;,;} in (5.16), computed with the B-spline 
basis shown in Fig. 5.34, is 

{{40.5, —112., —104., 80.1, —55.2, —36.8, 24.9, —18.1, -16.2, 

=179; 061,949, 215, 143,491, 138.179, 169, 

— 12.8, 6.89, 17.3, 21.7, 28.9, 40.4, 58.1, 83.1, 111., 133.}, 

{—157., -152., —105., —82.2, —65., —52.9, —44.7, —39., —36.3, 

— 4.42, 30.3, 30.4, 30.3, 28.3, 22.1, 9.55, 0.148, —1.41, 

5.3, 21.3, 17., 14.2, 13.2, 12.8, 12., 9.07, —1.03, —33.6}}. 


5.5 Dirichlet Problem in an Ellipse: Computational Difficulties 


5.5.1 Summary 


The boundary value problem and its analytic setup are the same in this example as 
in Sect. 5.2. The difference is that here we investigate the effect of both extreme and 
minor ill-conditioning, and discuss potential procedures for mitigating it. 
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5.5.2 Problem Statement 


Domain boundary. The curve OS is the ellipse parameterized by 
x1] =1+Cos[at], x2[t] = 4Sin[ar]t+5, O<t<2. (5.17) 


Its graph can be seen in Fig. 5.45. 


Governing equations. The displacement wu is the solution of the boundary value 
problem 


Z,oulx] =0, xES, 
ulx] = AX], xe. 


Test solution. The two components of the test solution u given by (5.3) are graphed 
in Fig. 5.40. 
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Fig. 5.40 The components of u[x], x € S. 


Boundary data function. The function Y is computed from the test solution u as 


Px] =ulx], xe OS. 


The graphs of its two components are displayed in Figs. 5.41 and 5.42, in Cartesian 
coordinates and parametric form, respectively. 
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Fig. 5.41 The components of [x] (Cartesian coordinates). 
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Fig. 5.42 The components of Y[x[t]] (parametric form). 


5.5.3 Solution Procedure 


Method. We apply the direct method, which reduces the problem to the boundary 
integral equation 


Vo() = (Wo+51)P onds, 


coded as 
¢ Dkrloobld= f Puyo Abldh +4 Aly 6.18) 
TWeak(x] Topv is] 


This equation is solved numerically to obtain an approximate density @. In turn, @ 
is used in the representation 


us]= $DbxsJooblals— f Pirlo Ablah, xes 6.19) 
Tr Tr 


to generate an approximate solution @ in S. 


Density function. The exact density @ is the boundary stress vector Tu computed 
from the test solution u and given by (5.4); its components are graphed in Figs. 5.43 
and 5.44. 


Fig. 5.43 The components of @[x] (Cartesian coordinates). 
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Fig. 5.44 The components of ¢[x[t]] (parametric form). 


Numerical approximation. We compute @ by the collocation method with a 
B-spline basis of elements b;,;. Then the approximate density is sought in the form 


LL sbi 
erat. | td 
l[x(¢]] = Ld ei, ;51,;[t] 


ij 


(5.20) 


where the numerical coefficients cg,;,; are determined by substituting (5.20) in 
(5.18). 

Since @ is continuous on the boundary, we use a piecewise linear spline with 
continuity at all knot locations. In terms of the parameterization (5.17), the knots 
and the smoothness at their locations are specified as the set 


1.1 3) -1 53-7 9 5 ll 3 13 7 15 
{0,848 398 1,354) 5) Bs a) B12h 
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Fig. 5.45 Left: the knot locations on 0S. Right: the collocation points. 


The 16 knots determine 16 functions b; ; for each component of @, generating a 
total of 32 basis functions. The graphs of the b;,; in parametric form are shown in 
Fig. 5.46. 

To compute @, we must choose an appropriate number of collocation points, 
which can be placed anywhere because D, P, @, and x are smooth on OS. The set of 
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16 values of the parameter f at our selection is 
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These values are marked on the graph on the right in Fig. 5.45. 
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Fig. 5.46 The 16 B-spline basis functions ); ;[f], O<t <2. 
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Conveniently, the collocation points are placed midway between the spline knots, 
where the 5; ; are not just continuous, but analytic. We use the same 16 points for 
both components of @, yielding a system of 32 constraining equations that enables 
us to compute the coefficients cg, ;,; in (5.20). The 32 x 32 matrix of this system has 
a condition number of 28,198. 


5.5.4 Solution 


Approximate density. The two components of @ are graphed in parametric form 
in Fig. 5.47, where the spline knot locations are also shown. 
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Fig. 5.47 The components of @[x[t]] (parametric form) and the knot locations. 


Ill-conditioning. The graphs in Fig. 5.47 bear no resemblance to those of the exact 
density ¢ displayed in Fig. 5.44. The first indication that something was wrong was 
the high condition number 28,198 of the coefficient matrix, which predicted extreme 
ill-conditioning. Also, the graphs themselves show classic signs of ill-conditioning. 
The singular value decomposition of the coefficient matrix contains two values that 
are essentially zero. As a consequence, the rank of the matrix is 2 units less than full 
rank, so the matrix is singular. The problem is caused by a combination of extreme 
symmetry in the domain, basis functions, and collocation points. This symmetry is 
broken if we move the collocation points away from the midpoint between the spline 
knots. 

We select new collocation points by shifting the previous ones by a very small 
amount (in this case, 1/1000) away from their original positions. Obviously, the 
new set 


127. 377) «=627)s 877_—s:«1127'- 1377: : 1627-1877 
2000? 2000° 2000° 2000° 2000’ 20007 20007 2000” 


2127 2377 2627 2877 3127 3377 3627 3877 
2000’ 2000’ 2000’ 2000’ 2000’ 2000’ 2000’ 2000 
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looks to the naked eye essentially the same as the original one on the right in Fig. 
5.45, but the modified numerical problem has a coefficient matrix with a condi- 
tion number of 673. Although the components of the new approximate density @, 
graphed in Fig. 5.48, still show signs of ill-conditioning, the true solution is begin- 
ning to emerge. 
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Fig. 5.48 The components of the new @[x{f]], (parametric form) and the knot locations. 


In Sect. 5.2, the collocation points for this problem were chosen to coincide with 
the spline knots, as opposed to the midpoint locations. That choice generated a coef- 
ficient matrix with a condition number of 24, meaning that the ill-conditioning had 
been removed and that the corresponding solution was very accurate. 

In Sect. 5.3, the same problem was solved with piecewise constant splines. 
There, the collocation points were positioned midway between the spline knot loca- 
tions, which resulted in good conditioning and a good outcome. For that example, 
choosing the collocation points to be the spline knots would produce the same ill- 
conditioning as exhibited in the piecewise linear spline version with midpoint loca- 
tions. The appropriate choice for piecewise constant splines is exactly the opposite 
of the appropriate choice for piecewise linear splines. This shows that there are no 
simple rules applicable in all circumstances. 

When ill-conditioning is suspected, another approach to deal with the situation is 
to create an overdetermined system by adding more collocation points. To illustrate 
this, we double the number of collocation points in the previous case from 16 to 32, 
choosing the new 32 points to correspond to the set of ¢ values 


{x 27 79 13 129 77 179 51 229 127 279 19 329 177 379 
100° 40? 200° 400° 50° 400° 200° 400° 100° 400° 200° 400° 25° 400° 200° 400° 


101 429 227 479 63 529 277 579 151 629 327 679 44 729 377 779 
100’ 400° 200’ 400’ 50’ 400° 200’ 400’ 100’ 400’ 200’ 400° 25° 400° 200° 400 


Keeping the same B-spline basis and knots but using the new set of collocation 
points yields an overdetermined system with a 64 x 32 coefficient matrix that has 
a condition number of 22. The graphs in Fig. 5.49 represent the new approximate 
density @ constructed from the solution of this system with the pseudo-inverse, and 
include the spline knot locations. 
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Fig. 5.49 The components of the overdetermined @|x{t]] (parametric form) and the knot locations. 
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Error analysis. This refined approximation @ agrees very well with the exact den- 
sity @ shown in Fig. 5.48. To visualize the difference between the two, in Fig. 5.50 
we graphed @ — @ relative to the maximum absolute value of @ on 0S. The same 
was done in parametric form in Fig. 5.51, which shows that the relative error in the 
two components of @ is only about 1%, except possibly on the extreme right side of 
OS, where the relative error can exceed 5%. 
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Fig. 5.51 The components of the relative error in @[x[t]] (parametric form). 
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5.6 Remarks. (i) The aim of the above computations is to show that ill-conditioned 
problems can arise quite naturally from seemingly appropriate choices. It would 
be logical, when using piecewise continuous linear spines, to choose the colloca- 
tion points away from the spline knots, where the basis functions do not have a 
continuous first-order derivative. But this choice is the very one that creates the ill- 
conditioning shown here. Also, it is quite possible to make choices that produce 
only moderate ill-conditioning. This situation is the most difficult to deal with as 
moderate ill-conditioning is often hard to separate from normal computational inac- 
curacy. Therefore, it is strongly recommended that careful analysis and vigilance be 
a part of every numerical problem. 

(ii) The set of coefficients {c1 jj, c2;,;} in (5.20), computed with the B-spline 
basis shown in Fig. 5.46 for the refined version of @[x{f]], is 


{{41.3, —69.9, -89.7, 73.7, —51.5, —34.4, —23.2, 
— 15.5, —4.49, 8.82, 19., 32., 51.4, 80., 111., 124.}, 
{—178., -191., —144., —100., —65.5, —40.5, —21.7, 
—4.51, 15.1, 22.3, 20.9, 19.3, 16.7, 8.91, —13.9, —74.}}. 


5.6 Dirichlet Problem in an Ellipse: Error Evaluation 


5.6.1 Summary 


We revisit the example in Sect. 5.2 using a piecewise linear spline with equally 
spaced knots and collocation points, but this time place the emphasis on evaluating 
the size of the error. 


5.6.2 Problem Statement 


Domain boundary. The curve OS is the ellipse parameterized by 
x1[t]=14Cos[at], x2[f] =5Sin[re]+3, O<1<2. (5.21) 


Its graph can be seen in Fig. 5.57. 


Governing equations. The displacement wu is the solution of the boundary value 
problem 


Zyoul4] =0, x ES, 
ulx]}= Ylx], xe. 
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Test solution. The two components of the test solution u given by (5.3) are graphed 
in Fig. 5.52. 
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2 
Fig. 5.52 The components of u{x], x € S. 
Boundary data function. The function Y is computed from the test solution ux] 
as 
Px] =ulx], xe OS. 


The graphs of its two components are displayed in Figs. 5.53 and 5.54, in Cartesian 
coordinates and parametric form, respectively. 


Fig. 5.53 The components of [x] (Cartesian coordinates). 
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Fig. 5.54 The components of Y[x[r]] (parametric form). 
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5.6.3 Solution Procedure 


Method. We apply the direct method, which reduces the problem to the boundary 
integral equation 


Vo(~) =(Wo+5 1A ond, 


coded as 
¢ Desloobldy= fp PleyoPbhlah+}P bi. (5.22) 
TWeak{x] TCPV x] 


This equation is solved numerically to obtain an approximate density @. In turn, @ 
is used in the representation 


ul] = $ Di ylo@hldty— $ PxyloPbldh, xs (5.23) 
Tr Tr 


to generate an approximate solution @ in S. 


Density function. The exact density @ is the boundary stress vector Tu computed 
from the test solution u and given by (5.4); its components are graphed in Figs. 5.55 
and 5.56. 
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Fig. 5.56 The components of ¢[x[t]] (parametric form). 
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Numerical approximation. We compute an approximate density @ by the collo- 
cation method with a B-spline basis of elements b;,;. Then the approximate density 
is sought in the form 


Lr iii] 


[x(t] = TTer4 jbiylf : (5.24) 
ij 


where the numerical coefficients cq ;,; are determined by substituting (5.24) in 
(5.22). 

Since @ is continuous on the boundary, we use a piecewise linear spline with 
continuity at all knot locations. In terms of the parameterization (5.21), the knots 
and the smoothness at their locations are specified as the set 


123 4 678 9 
{0, §) $39 1, $3 $9 5) 392} 
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Fig. 5.57 Left: the knot locations on 0S. Right: the collocation points. 


The 10 knots determine 10 functions b;,; for each component of @, generating a 
total of 20 basis functions. The graphs of the b;,; in parametric form are shown in 
Fig. 5.58. 

To compute ©, we must choose an appropriate number of collocation points, 
which can be placed anywhere because D, P, ~, and x are smooth on 0S. We choose 
these points to coincide with the knots, so the set of values of the parameter ¢ at our 
selection is 


These values are marked on the graph on the right in Fig. 5.57. 

The same 10 points are used for both components of @, yielding a system of 20 
constraining equations that enables us to compute the coefficients cq j,; in (5.24). 
The 20 x 20 matrix of this system has a condition number of 9. 


5.6 Dirichlet 


Problem in an Ellipse: Error Evaluation 
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Fig. 5.58 The 10 B-spline basis functions b,, ;[t], 0<t < 2. 
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Approximate density. The two components of @ are graphed in Figs. 5.59 and 


5.60. 


2 


Fig. 5.59 The components of @[x] (Cartesian coordinates). 
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Fig. 5.60 The components of @[x[t]] (parametric form) and the knot locations. 


We can now compare @ with the exact density @ displayed in Figs. 5.4 and 5.5. 


Error analysis. To visualize the difference between the approximate and exact den- 
sities, in Fig. 5.61 we graphed the function @ — @ relative to the maximum absolute 
value of @ on 0S. The same was done in parametric form in Fig. 5.62, which in- 
dicates that the relative error in the two components of @ is only about 5%, except 
possibly at the extreme right side of the boundary, where it can exceed 15%. 
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Fig. 5.62 The components of the relative error in @[x{¢]] (parametric form). 
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5.6.5 O(h’) Analysis 


It is well known that the optimal order of approximation for piecewise linear splines 
used in interpolation is O(h7), where h is the spacing between adjacent knots. Here, 
we conduct a computational study to measure the order of convergence in the di- 
rect boundary element method for the Dirichlet problem, which reduces to equation 
Di . The procedure consists in analyzing the rate of improvement of the error as h 
decreases. Since a global measurement of the error is required, we perform all the 
error measurements in the L?-norm. For the example in this section, using our 10 
basis functions with h = 0.2, we obtain the component error values 


|G: — Pi ||2 = 9.74, 
|| 2 — P2||2 = 6.14. 


We did a sequence of computational runs with an increasing number of basis func- 
tions, namely 


20, 30, 40, 50, 60, 70, 80, 90, 100, 


and calculated the L?-norm error for the two components in each case. The results 
are displayed in Fig. 5.63. 
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Fig. 5.63 The values of || @ — @q||2 as functions of basis size. 


As expected, the error decreases when the number of basis functions increases, 
which corresponds to a decrease in the size of h. To see if this decrease is O(h7), we 
attempt to fit the model 


| — pll2 = Ch" 
to the error. Taking logarithms on both sides yields 


Log|||@ — g||2] = Log|C] + nLog[h}. 
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Figure 5.64 shows the data points for Log|||@q — @q||2] as functions of Log[A]. 


Fig. 5.64 The values of Log|||@a — @a||2] as functions of Log[h]. 


The least-squares fit to this data leads to the values 


Log||| $1 — 91 ||2] & 5.25 + 2.12Log[h], 
Log|||@2 — @2||2] = 5.09 + 2.11Log{h]. 


When these models, superimposed on the data points in Fig. 5.64, are expressed in 
traditional form, they produce the approximations 


\| G1 — gi|l2 = 196n7!?, 
|| 2 — @2||2 = 166n7"". 


Our analysis strongly supports the conjecture that the collocation method with 
piecewise linear splines is indeed a second-order method. 


5.7 Remarks. (i) The above computational experiment has been performed for only 
one boundary integral method (the direct method), one boundary element method 
(the collocation method), and one boundary value problem (the Dirichlet problem). 
Additional work is needed to see if the result can be generalized to other choices of 
potential applications. 


(ii) The set of coefficients {cj j,;, c2,;,;} in (5.24), computed with the the B-spline 
basis shown in Fig. 5.58, is 


{{42.3, 96.2, —69.4, —36.7, —20., —4.57, 14.9, 34.4, 72.5, 129.}, 
{—179., —168., —90.8, —44.2, —14.8, 15.5, 22.4, 19., 12.8, —28.7}}. 
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5.7 Dirichlet Problem in a Domain with Corners: Piecewise 
Cubic Spline 


5.7.1 Summary 


This example illustrates the application of the direct method to a Dirichlet problem 
in a domain with a continuous and asymmetric boundary that has three corners. The 
numerical procedure is, in general, the same as in the example discussed in Sect. 
5.2, but here the method is implemented by means of a piecewise cubic spline with 
equally spaced knots and collocation points, and basis functions with discontinuities 
at the knots. 


5.7.2 Problem Statement 


Domain boundary. QS is the curve parameterized by 


1 (3-+5¢+487? —32), O<7¢< 1, 


12 
x1 [t] = ¢ 5 (15-78), Le¢e2., 
4 (5-21), Perey, 
(5.25) 
at, 0<7< 1, 
x2[t] = ¢ ge (86 — 1274 + 1282? — 3277), 1<t<2, 
3-1, 2<t<3. 


Its graph can be seen in Fig. 5.70. 


Governing equations. The displacement u is the solution of the boundary value 
problem 


Zyoulx] =0, xES, 
ulx]}= Pix], xed. 


Test solution. The two components of the test solution u given by (5.3) are graphed 
in Fig. 5.65. 


Boundary data function. The function Y is computed from the test solution wu as 
P\x| =ulx], xe Os. 


The graphs of its two components are displayed in Figs. 5.66 and 5.67, in Cartesian 
coordinates and parametric form, respectively. 
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Fig. 5.65 The components of u|x], x € S. 
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Fig. 5.66 The components of Y |x] (Cartesian coordinates). 
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Fig. 5.67 The components of [x(t]] (parametric form). 
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5.7.3 Solution Procedure 


Method. We apply the direct method, which reduces the problem to the boundary 
integral equation 


Vo(~) =(Wo+5 1A ond, 


coded as 
¢ Desloobldy= fp PleyoPbhlah+}P bi. (5.26) 
TWeak{x] TCPV x] 


This equation is solved numerically to obtain an approximate density @. In turn, @ 
is used in the representation 


ul] = $ Di ylo@hldty— $ PxyloPbldh, xs 65.27) 
Tr Tr 


to generate an approximate solution @ in S. 


Density function. The exact density @ is the boundary stress vector Tu computed 


from the test solution wu and given by (5.4); its components are graphed in Figs. 5.68 
and 5.69. 


Fig. 5.69 The components of ¢[x[t]] (parametric form). 
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Numerical approximation. We compute @ by the collocation method with a 
B-spline basis of elements b;,;. Then the approximate density is sought in the form 


Ldc1 i, /bi,;{t] 
x 


@[x{t]] = Vdc ashi (¢] ’ (5.28) 
J 


1 


where the numerical coefficients cq; are determined by substituting (5.28) in 
(5.26). 

Since @ is discontinuous at the three corners of the boundary, the b; ; need to 
exhibit the same feature. In our case, we choose a piecewise cubic spline that is 
twice continuously differentiable at the secondary knots and discontinuous at the 
primary knot locations t = 0, 1, 2,3. In terms of the parameterization (5.25), the 
primary and secondary knots and the smoothness at their locations are specified as 
the set 


{0,0,0,04,1, 1,1, 1, 3, 2,2, 2, 2, 3, 3, 3,3, 3}. 


The knots are marked on the graph on the left in Fig. 5.70. 
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Fig. 5.70 Left: the knot locations on 0S. Right: the collocation points. 


The 6 knots determine 15 functions b;,; for each component of @, generating a 
total of 30 basis functions. The graphs of the b;,; in parametric form are shown in 
Fig. 5.71. 

To compute @, we must choose an appropriate number of collocation points, 
which can be placed anywhere on 0S except at the corners t = 0, 1, 2,3, where @ is 
expected to be discontinuous. The set of 15 values of the parameter f at our selection 
is 


1 317 9 Ht 13 3 17 19 21 23 5 27 29 
10° 10° 2° 107 10° 10° 10’ 2° 10’ 10° 10° 10° 27 10° 10 S° 


These values are marked on the graph on the right in Fig. 5.70. 

The same 15 points are used for both components of @, yielding a system of 30 
constraining equations that enables us to compute the coefficients cq j,; in (5.28). 
The 30 x 30 matrix of this system has a condition number of 124. 
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Fig. 5.71 The 15 B-spline basis functions b,, ;[t], 0<t < 3. 
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5.7.4 Solution 


Approximate density. The two components of @ are graphed in Figs. 5.72 and 
5.73. The latter also shows the spline knot locations. It can be seen that @ agrees 
very well with the exact density @ = Tu in Figs. 5.68 and 5.69. 


—100 - 


Fig. 5.73 The components of @[x[t]] (parametric form) and the knot locations. 


Approximate solution. We use @ and # in (5.27) to compute an approximation 
ii to the exact solution u in S. The graphs of the two components of @# are shown in 


Fig. 5.74. 


N 


2 


Fig. 5.74 The components of di[x], x € S. 
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Error analysis. To visualize the difference between the approximate and exact 
densities, in Fig. 5.75 we graphed @ — @ relative to the maximum absolute value of 
@ on OS. The same was done in parametric form in Fig. 5.76, which indicates that 


the relative error in the two components of @ is only about 1%, except possibly at 
the extreme right segment of the boundary, where it can exceed 5%. 
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Fig. 5.76 The components of the relative error in @[x[¢]] (parametric form). 


5.8 Remark. The set of coefficients {c1 jj, c2,;,;} in (5.28), computed with the B- 
spline basis shown in Fig. 5.71, is 


{{13.4, 25.6, 33.4, 149., 169., —157., —98.6, —41.2, 
25.5, -18.5, —22.8, 


16.6, —4.93, 9.33, 15.}, 
{-10.9, 6.8, 6.86, 45.8, —44.7, —133., —114., 


72.1, 49.6, —39.3, —14., —4.22, 19.5, 21.4, 20.1} }. 
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5.8 Dirichlet Problem in a Domain with Corners: Piecewise 
Quintic Spline 


5.8.1 Summary 


The analytic setup of the problem in this example is the same as in Sect. 5.7, with 
the same asymmetric boundary curve 0S. We show how the accuracy of the method 
can be improved even when the boundary has corners if we increase the degree of 
the spline as well as the number of collocation points. We implement the colloca- 
tion method using a piecewise fifth-degree polynomial spline over equally spaced 
knots and collocation points, and discuss the accuracy and ill-conditioning of the 
procedure. 


5.8.2 Problem Statement 


Domain boundary. QS is the curve parameterized by 


2 (3+5¢+48r? —3277), O<¢<1, 


12 
x1 [t] = ¢ 5 (15-72), 1<¢< 2,5 
4 (5-21), 22S, 
(5.29) 
at, he aoe 0 
x2[t] = ¢ ge (86 — 1274 + 1282? — 3273), 1<t<2, 
3-1, 2<t<3. 


Its graph can be seen in Fig. 5.82. 


Governing equations. The displacement u is the solution of the boundary value 
problem 


Zyoulx] =0, xES, 
ulx]}= Pix], xed. 


Test solution. The two components of the test solution u given by (5.3) are graphed 
in Fig. 5.77. 


Boundary data function. The function Y is computed from the test solution u as 
P\x| =ulx], xeos. 


The graphs of its two components are displayed in Figs. 5.78 and 5.79, in Cartesian 
coordinates and parametric form, respectively. 
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Fig. 5.78 The components of Y [x] (Cartesian coordinates). 
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Fig. 5.79 The components of Y[x[t]] (parametric form). 
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5.8.3 Solution Procedure 


Method. We apply the direct method, which reduces the problem to the boundary 
integral equation 
Vo(p) = (Wo+51)P onds, 


coded as 


¢ Dkslooblaz= f Peso Pbldhy +4 Al. 


TWeak(x] 


(5.30) 
Topvpy 


This equation is solved numerically to obtain an approximate density @. In turn, @ 
is used in the representation 


ul] = $ Di slo@hldty— $ Px yloPbldh, xes 631) 
F Tr 


to generate an approximate solution @ in S. 
Density function. The exact density @ is the boundary stress vector Tu computed 


from the test solution u and given by (5.4); its components are graphed in Figs. 5.80 
and 5.81. 


Fig. 5.81 The components of ¢[x[t]] (parametric form). 
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Numerical approximation. We compute @ by the collocation method with a 
B-spline basis of elements b;,;. Then the approximate density is sought in the form 


Ldc1 i, ji, {t] 
zz 


@[x{t]] — Yeoh (t] ’ (5.32) 
ij 

where the numerical coefficients cg j;,; are determined by substituting (5.32) in 

(5.30). 

Since @ is discontinuous at the three corners of the boundary, the b; ; need to 
exhibit the same feature. In our case, we choose a piecewise quintic spline that is 
four times continuously differentiable at the secondary knots and discontinuous at 
the primary knot locations t = 0, 1, 2, 3. In terms of the parameterization (5.29), the 
primary and secondary knots and the smoothness at their locations are specified as 
the set 

{0, 0, 0,0, 0,0, 5, 4,3, 7,3, $,1,1,1,1,1,1, 8, 3, 2, 


11 12 13 15 16 17 18 19 20 
G2, 3,2,2,2,2,2,2, 3,8, 7,3, oO 3,3,3,3,3, 3}. 


The knots are marked on the graph on the left in Fig. 5.82. 


Fig. 5.82 Left: the knot locations on 0S. Right: the collocation points. 


The 22 knots determine 36 functions b;,; for each component of @, generating a 
total of 72 basis functions. The graphs of the b;,; in parametric form are shown in 
Fig. 5.83. 

To compute @, we must choose an appropriate number of collocation points, 
which can be placed anywhere except at the corners t = 0, 1, 2, 3, where @ is ex- 
pected to be discontinuous. The set of 72 values of the parameter f at our selection 
is 


{2 1 5 7 3 Wt 13 5 17 19 7 23 25 9 29 31 Il 35 37 
48° 16? 48° 48° 16? 48? 48° 16? 48? 48° 16? 48° 48° 16? 48° 48° 16? 48° 48° 
13 41 43 15 47 49 17 53 55 19 59 61 21 65 67 23 71 73 25 
16’ 48° 48? 16’ 48° 48? 16? 48° 48? 16? 48° 487 16? 48° 48? 16? 48? 48° 16? 
71 79 27 83 85 29 89 91 31 95 97 33 101 103 35 107 109 37 
48° 48> 16? 487 48> 16? 487 48° 16? 48° 48° 16? 48° 48 16? 48° 48° 16° 
113 115 39 19 121 41 y95q4g 127 43 131 133 45 137 139 47 143 
48° 48 16? 48° 48> 16? 1 48 16> “48> 48> 16> 48 48> 167 48 


These values are marked on the graph on the right in Fig. 5.82. 
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Fig. 5.83 The 36 B-spline basis functions b,, ;[t], 0<t < 3. 
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The same 72 points are used for both components of @, yielding a system of 144 
constraining equations that enable us to compute the coefficients cg; ; in (5.32). The 
144 x 72 matrix of this overdetermined system has a condition number of 712. 


5.8.4 Solution 


Approximate density. The two components of @ are graphed in Figs. 5.84 and 
5.85. The latter also shows the spline knot locations. It can be seen that @ agrees 
very well with the exact density @ = Tu in Figs. 5.80 and 5.81. 
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Fig. 5.85 The components of @[x[t]] (parametric form) and the knot locations. 


Approximate solution. We use @ and # in (5.31) to compute an approximation 
i to the exact solution u in S. The graphs of the two components of @ are shown 
in Fig. 5.86. Of course, to the naked eye they look no different than the graphs 
of the exact solution in Fig.5.74, which were drawn based on computation with a 
piecewise cubic spline. 
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Fig. 5.86 The components of ix], x € S. 


Error analysis. To visualize the difference between the approximate and exact 
densities, in Fig. 5.87 we graphed @ — @ relative to the maximum absolute value 
of @ on OS. The same was done in parametric form in Fig. 5.88, which indicates 
that the relative error in the two components of @ is only about 0.005%, except, as 
expected, at the corners of the boundary. 


Fig. 5.87 The components of the relative error in @[x] (Cartesian coordinates). 
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Fig. 5.88 The components of the relative error in @[x{t]] (parametric form). 
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This high level of accuracy illustrates the benefit of using a piecewise higher- 
degree polynomial spline for problems with smooth boundary conditions between 
adjacent corners. The accuracy in the vicinity of the corners is expected to deteri- 
orate significantly. We can partially mitigate this effect by using unequally spaced 
spline knots and carefully chosen collocation points, as we have done here. 


Ill-conditioning. We carried out a second numerical experiment (not shown here) 
with the collocation points placed midway between the spline knots. This produced 
a highly ill-conditioned problem. Consequently, we followed the procedure outlined 
in Sect. 5.5 and augmented the number of collocation points, which reduced from 
five to three digits the condition number for the coefficient matrix of the system that 
determines the cq,;,;, and thus led to an acceptable computational scheme. 


5.9 Remark. The set of the coefficients {cj jj, c2,i,;} in the approximation @ given 
by (5.32) for the basis functions shown in Fig. 5.83 is 
{{10.8, 15.7, 19.7, 23.8, 31., 46.6, 71.8, 108., 142., 
163., 160., 136., —157., —146., —128., —104., —80., —58.6, 
43.3, —32.4, —26., —22.1, —19.8, —18.7, —23.2, —21.9, —19.5, 
16.3, —12.4, —7.23, —1.43, 4.32, 8.67, 12., 14.3, 15.6}, 
{—13.2, —7.99, —2.01, 2.54, 6.98, 13.1, 20.3, 27.5, 26.2, 
8.12, —38.4, —87.3, —134., —131., —124., —113., —99.2, —83.5, 
69.9, —58.2, —50.2, —44.6, —41.1, —39.4, —14.5, —12.5, —8.71, 
— 3.28, 3.58, 11.4, 17.2, 20., 20.6, 20.7, 20.6, 20.6} }. 
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5.9.1 Summary 


This example illustrates the application of the boundary element method in a case 
where the boundary conditions are not constructed from a test solution. Conse- 
quently, the final result of our computation is not known to us a priori, which raises 
the question of validation. We aim to answer this question in a ‘natural’ way. Here, 
the domain S is a square with Dirichlet conditions prescribed on its boundary. We 
obtain an approximate solution by means of the direct method and a piecewise cubic 
spline on equally spaced knots and collocation points. In the next section, we use 
this solution as input for a Neumann boundary value problem, whose output will 
enable us to settle the validation issue for the current example. 
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5.9.2 Problem Statement 


Domain boundary. The contour 0S is the square parameterized by 


t, 0<t<1, 0, 0<t<1, 
1 1<t<2 t—-1l, 1<t<2 
xIfJ= eo” SES otf a (5.33) 
3-t, 2<t<3, 1, 2<t<3, 
0, 3<t<4, 4-t 3<t<4. 


Its graph can be seen in Fig. 5.91. 


Governing equations. The displacement u is the solution of the boundary value 
problem 


Zyoulx] =0, xES, 
ulx]}= Ylx], xe. 
Test solution. No test solution is used in this example. 


Boundary data function. The function Y is prescribed without reference or con- 
nection to any known test solution. Our choice, in parametric form, is 


5+ 1017, 0<r<1, 
Arpt) = 20—10t+5t7, 1<t<2, 
=124-32-8F,. 2<7<3, 
—$1+4-7F, G<t<4, 
15=10r; O=7< iy 
Parle] = 17-4 1, Tere, 
37 20r+-5r",, 2<2< 4, 


—41+42-— 7, 3< tf <4, 


The graphs of its two components are displayed in Figs. 5.89 and 5.90, in Cartesian 
coordinates and parametric form, respectively. 
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Fig. 5.89 The components of [x] (Cartesian coordinates). 
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Fig. 5.90 The components of Y[x[t]] (parametric form). 


5.9.3 Solution Procedure 


Method. We apply the direct method, which reduces the problem to the boundary 
integral equation 


Vo(~) =(Wo+51)A ond, 
coded as 
Dinyo@pldh= f Plexo Pbhlah+}P7 bi. (5.34) 
TWeak(x| Topvpy 


This equation is solved numerically to obtain an approximate density @. In turn, @ 
is used in the representation 


ub] = PD sJooblaly— J Pirlo Pblah, xes 65.35) 
r E 


to generate an approximate solution @ in S. 


Density function. This is not known beforehand since we are not using a test solu- 
tion. 


Numerical approximation. We compute @ by the collocation method with a 
B-spline basis of elements b; ;. Then the approximate density is sought in the form 


Lhe, (i,j) 
~ ij 


@|x[t]] = F300, ;bi lt] (5.36) 
ij 
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where the numerical coefficients cq ;; are determined by substituting (5.36) in 
(5.34). 

Since @ is expected to be discontinuous at the four corners of 0S, the bi, j need 
to exhibit the same feature. In our case, we choose a piecewise cubic spline that is 
twice continuously differentiable at the secondary knots and discontinuous at the 
primary knot locations t = 0, 1, 2, 3, 4. In terms of the parameterization (5.33), the 
primary and secondary knots and smoothness at their locations are specified as the 
set 


{0, 0, 0,0, 5,1, 1,1, 3, 2,2, 2, 3,3, 3,3, 4,4,4,4, 4}. 


The knots are marked on the graph on the left in Fig. 5.91. 


Fig. 5.91 Left: the knot locations on 0S. Right: the collocation points. 


The 9 knots determine 20 functions b; ; for each component of @, generating a 
total of 40 basis functions. The graphs of the b;,; in parametric form are shown in 
Fig. 5.92. 

To compute @, we must choose an appropriate number of collocation points, 
which can be placed anywhere except at the corners t = 0, 1, 2,3,4, where ¢— is 
expected to be discontinuous. The set of 20 values of the parameter f at our selection 
is 

1 3 17 9 UW 13 3 17 19 21 23 5 27 29 31 33 7 37 39 
10’ 10° 2° 10° 10° 10° 10° 2° 10° 10° 10° 10° 2° 10° 10° 10> 10° 2? 10> 10S" 


These values are marked on the graph on the right in Fig. 5.91. 

The same 20 points are used for both components of @, yielding a system of 40 
constraining equations that enables us to compute the coefficients cq j,; in (5.36). 
The 40 x 40 matrix of this system has a condition number of 120. 
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Fig. 5.92 The 20 B-spline basis functions b,, ;[t], O<t <4. 
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5.9.4 Solution 


Approximate density. The components of @ are graphed in Figs. 5.93 and 5.94. 


1 


Fig. 5.93 The components of @[x] (Cartesian coordinates). 


Fig. 5.94 The components of @[x[t]] (parametric form). 


Approximate solution. We use @ and # in (5.35) to compute an approximation 


ii to the exact solution u in S. The graphs of the components of # are shown in Fig. 
aS. 


Error analysis. The exact solution of this problem is not known; therefore, no 
comparison can be made between the calculated density @ and the exact density @. 
However, since @ represents Tu, in Sect. 5.10 we use the approximate density @ 
computed here as the input function to a Neumann boundary value problem; that is, 


2lx({¢]] = Tu) [ale]] = Glee] on a. 


The unknown density y in the integral equation for that problem represents the 
value of the solution u on 0S; hence, it is the (approximate) data function for the 
example discussed above, and we can compare it to our choice of function Y for 
validation purposes. 
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Fig. 5.95 The components of di[x], x € S. 


5.10 Remarks. (i) As already mentioned, @ is not constructed from a known test so- 
lution like (5.3). The family of test solutions that we generated in Sect. 4.13, which 
includes (5.3), exhibits radial symmetry about the point x in polar coordinates, so 
it is a very limited collection of solutions for the general problem. The paramet- 
ric functions used to define Y in this example does not have such symmetry and, 
therefore, it is a good illustration of the general case. 


(ii) The set of coefficients {c1 jj, c2;,;} in (5.36), computed with the B-spline 
basis shown in Fig. 5.92, is 
{{—22.8, 18.2, —13.3, 12.6, 37.4, 91.3, 61.5, 44.5, 35.5, 23.4, 
18.9, 24.8, 28., 5.41, —18.7, 74.1, —36.9, -13., —3.61, —2.35}, 
{-55.6, -15:, 1.48, 6:75, 13,8, 369, 39.4, 96:7, 1.48, 19.6, 
117., 84.7, 54.1, 33.1, 13., 17.8, 25.4, 29.9, 6.16, —23.5}}. 


5.10 Neumann Problem in a Square: Validation Procedure 
5.10.1 Summary 


This example illustrates the application of the boundary element method to a Neu- 
mann problem in a square, where the function prescribed on the boundary is the 
approximate density computed for the Dirichlet problem in Sect. 5.9. The solu- 
tion is obtained with the direct method by means of a piecewise cubic spline on 
equally spaced collocation points. The numerical output will be compared to the 
input Dirichlet data in the preceding example, to help illustrate the validity of our 
methods when the exact solution is unknown. 
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5.10.2 Problem Statement 


Domain boundary. 0S is the square parameterized by 


ig 0<r<il, 0, 0<t<1, 
1 1<t<2 t-1, 1<t<2 
xiao Sheek si = (5.37) 
3-t, 2<t<3, 1, 2<t<3, 
0, 3<t<4, 4-t 3<t<4. 


Its graph can be seen in Fig. 5.98. 


Governing equations. The displacement u is the solution of the boundary value 
problem 


Zyoulx] =0, xES, 
(Tu) [x] = T.oulx] = Qix], xe OS. 


Test solution. No test solution is used for this example. 
Boundary data function. The function 2 is defined by 


Qhx]=Glx], O<t<4, 


where @ is the approximate density computed in Sect. 5.9. The graphs of the two 
components of 2 are displayed in Figs. 5.96 and 5.97, in Cartesian coordinates and 
parametric form, respectively. 

This choice makes it possible for us to verify the computational accuracy in the 
boundary element method used to calculate both @ and the approximate density w 
in this example. However, no assessment of the exact accuracy can be made because 
the exact solution is unknown. 


1 


Fig. 5.96 The components of 2[x] (Cartesian coordinates). 
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ee 


75 


Fig. 5.97 The components of 2[x/t]] (parametric form). 


5.10.3 Solution Procedure 


Method. We apply the direct method, which reduces the problem to the boundary 
integral equation 


(Wo + sl)w-VW2 on 0s, 


codes as 
¢ Prlovblay+svkl= ¢ DixxJo2plany. (5.38) 
TCpv ix] TWeak{x] 


This equation is solved numerically to obtain an approximate density . In turn, w 
is used in the representation 


ule] = § Disy]o2hldhy— fp Pe.rlowbldh, xes 65.39) 
Tr Tr 


to generate an approximate solution d in S. 


Density function. This example does not use a test solution u, so the exact density 
function y is not available. 


Numerical approximation. We compute by the collocation method with a 
B-spline basis of elements b;,;. Then the approximate density is sought in the form 


Lhe 4,j51,,[t) 
S ij 


Wixle]] = LDer, jbi,j{d] , (5.40) 
ij 


where the numerical coefficients cq ;,; are determined by substituting (5.40) in 
(5.38). 
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Since y represents u on OS and the latter is continuous but has discontinuous 
derivatives at the corners, the b;,; need to exhibit the same feature. In our case, 
we choose a piecewise cubic spline that is twice continuously differentiable at the 
secondary knots and continuous but with discontinuous derivatives at the primary 
knot locations t = 0, 1, 2, 3, 4. In terms of the parameterization (5.37), the knots and 
the smoothness at their locations are specified as the set 


{0, 0,0, 5,1, 1,1, 4, 2, 2,2, 3, 3, 3, 3, 4,4, 4,4}. 


The knots are marked on the graph on the left in Fig. 5.98. 


0 1 0 1 


Fig. 5.98 Left: the knot locations on 0S. Right: the collocation points. 


The 9 knots determine 16 functions b;,; for each component of W, generating a 
total of 32 basis functions. The graphs of the b;,; in parametric form are shown in 
Fig. 5.99. 

To compute W, we must choose an appropriate number of collocation points, 
which can be placed anywhere on 0S because y is expected to be continuous on the 
boundary. The set of 16 values of the parameter ¢ at our selection is 


(3 13 15 17 19 21 23 25 27 29 31 
8° 8° 8° 8° 8? 8° 8? 8? 8° 8° 8B? 8° 8° 8°? 8° 8 
These values are marked on the graph on the right in Fig. 5.98. 

The same 16 collocation points are used for both components of W, yielding a 
system of 32 constraining equations that enable us to compute the coefficients cg ;j, j 
in (5.40). The 32 x 32 matrix of this system has a condition number of 256. 


5.10 Neumann Problem in a Square: Validation Procedure 
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Fig. 5.99 The 16 B-spline basis functions bj, ;[t], 0<t <4. 
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5.10.4 Solution 


Approximate density. The two components of W are graphed in Figs. 5.100 and 
5.101. 


1 1 


Fig. 5.100 The components of [x] (Cartesian coordinates). 


1 2 3 4 


Fig. 5.101 The components of [x[t]] (parametric form). 


Approximate solution. We use % and @ in (5.39) to compute an approximation 
ii to the exact solution u in S. The graphs of the two components of @# are shown in 
Fig. 5.102. 


Error analysis. In Remark 5.10(i), we mentioned that y represents the boundary 
values of u. Since those values were given by the function Y in the preceding 
example, and since & has been calculated here using an input constructed from Y, 
it is important for us to determine how close both sets of results are. 
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Fig. 5.102 The components of @i[x], x € S. 


To visualize the difference between W and Y, in Fig. 5.103 we graphed w— FY 
relative to the maximum absolute value of Y on OS. These graphs indicate that the 
relative error in the two components of W is only about 0.2%, except possibly at 
the corners of the boundary, where it can exceed 0.3%. The fact that the composite 
of both the Dirichlet and Neumann problems yields such a small error supports the 
validity of our calculations. 


—0.002 + 


Fig. 5.103 The components of the relative error in W[x[t]] (parametric form). 


5.11 Remarks. (i) The exact solution of this problem is not known; therefore, no 
comparison can be made between the calculated density ~ and the exact density yw. 


(ii) The set of coefficients {c jj, c2;,;} in (5.40), computed with the B-spline 
basis shown in Fig. 5.99, is 
{{5., 5.02, 6.67, 11.6, 15., 15., 15.8, 18.3, 19.9, 19.9, 18.6, 
14.6, 11.9, 11.9, 10.8, 7.33}, {14.9, 14.9, 13.3, 8.33, 5.01, 5.01, 
7.01, 13., 17., 16.9, 17.8, 20.3, 21.9, 21.9, 20.8, 17.3}}. 
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5.11 Neumann Problem in a Square: Piecewise Cubic Spline 


5.11.1 Summary 


This example illustrates the application of the direct method to a Neumann prob- 
lem in a square domain, where the boundary condition is computed from a known 
test solution. We use the collocation method with piecewise cubic splines on equally 
spaced knots and collocation points to construct an approximation of the original so- 
lution from the computed boundary data function. The results are validated through 
a comparison between the exact and approximate solutions. 


5.11.2 Problem Statement 


Domain boundary. The contour 0S is the square parameterized by 


t, 0<t<1, 0, 0<t<1, 
1 bere 2, t-1, 1<t<2 
xI[fJ= eo we sold bt (5.41) 
3-t, 2<1<3, 1, 2<t<3, 
0, 3<t<4, 4-t 3<t<4. 


Its graph can be seen in Fig. 5.109. 
Governing equations. The displacement wu is the solution of the boundary value 
problem 

Zyou[x]=0, xES, 

(Tu) [x] = T,oulxy] = Qh], xe dS. 


Test solution. The two components of the test solution u given by (5.3) are graphed 
in Fig. 5.104. 


Fig. 5.104 The components of u[x], x € S. 
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Boundary data function. The function 2 = Tu is given by (5.4). The graphs of 
. 5.105 and 5.106, in Cartesian coordinates 


its two components are displayed in Figs 
and parametric form, respectively. 


1 


1 


Fig. 5.105 The components of 2[x] (Cartesian coordinates). 


\ 
\ 


Fig. 5.106 The components of 2[x[t]] (parametric form). 


5.11.3 Solution Procedure 
Method. We apply the direct method, which reduces the problem to the boundary 
integral equation 

(Wo + 51) w-VW2 onds, 


coded as 
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¢ Prlovblay+3vkl= ¢ DixyJo2plary. (5.42) 
Topvp] TWeak|x] 
This equation is solved numerically to obtain an approximate density . In turn, W 
is used in the representation 


ule] = f Dixy]o2hlaly— f P.rlowbldh, xes (5.43) 
Tr Tr 


to generate an approximate solution @ in S. 


Density function. The exact density y is the restriction to 0S of the test solution 
u; its components are graphed in Figs. 5.107 and 5.108. 


1 1 


Fig. 5.107 The components of y[x] (Cartesian coordinates). 


10 


Fig. 5.108 The components of w[x[t]] (parametric form). 
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Numerical approximation. We compute by the collocation method with a 
B-spline basis of elements b; ;. Then the approximate density is sought in the form 


Lde1 ji, lt] 
ai ij 


W[x[t]] = EEen,jbi lt] } (5.44) 
ij 

where the numerical coefficients cq ;,; are determined by substituting (5.44) in 

(5.42). 

Since y represents u on OS and the latter is continuous but has discontinuous 
derivatives at the corners, the b;,; need to exhibit the same feature. In our case, 
we choose a piecewise cubic spline that is twice continuously differentiable at the 
secondary knots and continuous but with discontinuous derivatives at the primary 
knot locations t = 0, 1, 2, 3, 4. In terms of the parameterization (5.41), the knots and 
the smoothness at their locations are specified as the set 


{0, 0, 0, 5 1, 1, 1, 3,2, 2, 2, 3, 3,3, 3, 5,4, 4,4}. 


The knots are marked on the graph on the left in Fig. 5.109. 


| 
| | 
| | 


0 1 0 1 


Fig. 5.109 Left: the knot locations on 0S. Right: the collocation points. 


The 8 knots determine 16 functions b;,; for each component of W, generating a 
total of 32 basis functions. The graphs of the b;,; in parametric form, modified to 
enforce the required continuity at the knot locations, are shown in Fig. 5.110. 

To compute W, we must choose an appropriate number of collocation points, 
which can be placed anywhere because y is expected to be continuous on 0S. We 
select equally spaced points that avoid the primary knots. The set of 16 values of the 
parameter ¢ at our selection is 


{BEEP eee ee pre ee Fz 


87> 87 8° 8° 8? 87 8? 8? 8°? 8° 8° 8° 8° 8° 8° 8 


These values are marked on the graph on the right in Fig. 5.109. 
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Fig. 5.110 The 16 B-spline basis functions b;,;[f], 0<t < 4. 


The same 16 collocation points are used for both components of W, yielding a 
system of 32 constraining equations that enable us to compute the coefficients cg j, ; 
in (5.44). The 32 x 32 matrix of this system has a condition number of 256. 
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5.11.4 Solution 


Approximate density. The two components of W are graphed in Figs. 5.111 and 
5.112. The latter also shows the spline knot locations. It can be seen that W agrees 
very well with the exact density y in Figs. 5.107 and 5.108. 


eee 


1 


Fig. 5.111 The components of W[x] (Cartesian coordinates). 


10 


Fig. 5.112 The components of [x[t]] (parametric form). 


Approximate solution. We use % and 2 in (5.43) to compute an approximation 
ii to the exact solution u in S. The graphs of the two components of a are shown in 
Fig. 5.113. 
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Fig. 5.113 The components of di[x], x € S. 


Error analysis. To visualize the difference between the approximate and exact 
densities, in Fig. 5.114 we graphed w — y relative to the maximum absolute value 
of yon ds. 


Fig. 5.114 The components of the relative error in W[x] (Cartesian coordinates). 


The same was done in parametric form in Fig. 5.115, which also shows the pri- 
mary and secondary knots and indicates that the relative error in the two components 
of W is about 0.05%. 
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Fig. 5.115 The components of the relative error in W[x[t]] (parametric form). 


With such a small error on the boundary, we expect an equally good agreement 
in the interior of the domain. This is indeed confirmed by Fig. 5.116, which displays 
the graphs of the components of #—u in S. 


0.001 


-0.001 


Fig. 5.116 The components of a[x] — ux], x € S. 


5.12 Remark. The set of coefficients {cj jj, c2;,;} in (5.44), computed with the 
B-spline basis shown in Fig. 5.110, is 


{{0.474, 0.493, 0.523, 0.632, 0.706, 1.84, 4.4, 8.07, 

10.3, 8.9, 6.62, 5.15, 4.52, 3.45, 1.82, 0.864}, 

{—-14.8, 16.4, —20.8, —28.8, —34.2, —34.8, —36.7, —40.5, 
43., —36.4, —26.6, —21.1, —18.9, —17.7, —15.8, —15.}}. 
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5.12 Robin Problem in a Square: Piecewise Cubic Spline 


5.12.1 Summary 


Here, we illustrate the application of the direct method to the Robin boundary value 
problem in a square domain. We compute the boundary condition from a known test 
solution, then use the collocation method with a piecewise cubic spline on equally 
spaced knots and collocation points to approximate the exact solution. The results 
are validated through a comparison between the approximate and exact solutions. 


5.12.2 Problem Statement 


Domain boundary. The contour 0S is the square parameterized by 


t, 0<t<1, 0, 0<t<1, 
1 1<t<2 t-1, 1<t<2 
xI[fJ= SES gies aia (5.45) 
3-t, 2<1<3, 1, 2<t<3, 
0, 3<t<4, 4-t 3<t<4. 


Its graph can be seen in Fig. 5.122. 


Governing equations. The displacement u is the solution of the boundary value 
problem 


Zyoulx] =0, xES, 
(Tu) [x] + oul[x] = 4 |x], x€ OS. 


Test solution. The two components of the test solution u given by (5.3) are graphed 
in Fig. 5.117. 


Fig. 5.117 The components of u[x], x € S. 
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Boundary data function. The function .% is given by 
(Tu) [x]+ oul] = “[x], xe OS, 


where, for computational purposes, we have chosen the particular 2 x 2 matrix 


The graphs of the two components of .% are displayed in Figs. 5.118 and 5.119, in 
Cartesian coordinates and parametric form, respectively. 


1 1 


Fig. 5.118 The components of .% [x] (Cartesian coordinates). 


meena ll 
—150 


Fig. 5.119 The components of .% [x[t]] (parametric form). 
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5.12.3 Solution Procedure 


Method. We apply the direct method, which reduces the problem to the boundary 
integral equation 


(Wo+ 51 w+Vo(ow) =Vo(%) onds, 


coded as 


¢ Pslowblah+ ¢ Dkxylo(owb))ah +} vb 


TCpv ix] TWeak{x] 
a $ Dix, y]o # [y|dG. (5.46) 
TWeak{x] 


This equation is solved numerically to obtain an approximate density W. In turn, w 
is used in representation 


ul[x] = f Dixy 0 (4H — owly]) di, — f Phx,yI owb]dIy, xEs (5.47) 
r r 


to generate an approximate solution @ in S. 


Density function. The exact density y is the restriction to 0S of the test solution 
u; its components are graphed in Figs. 5.120 and 5.121. 


1 1 


Fig. 5.120 The components of wx], x € 0S (Cartesian coordinates). 
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10 
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Fig. 5.121 The components of y{x{t]], 0 < t <4 (parametric form). 


Numerical approximation. We compute by the collocation method with a 
B-spline basis of elements b;,;. Then the approximate density is sought in the form 
Then the approximate density is sought in the form 


rei jit] 
a ae | 


Wlxle]] _ Ldc25,;bi,{t] ? (5.48) 
ij 
where the coefficients cq ;,; are determined by substituting (5.48) in (5.46). 

Since y represents u on OS and the latter is continuous but has discontinuous 
derivatives at the corners, the b;,; need to exhibit the same feature. In our case, 
we choose a piecewise cubic spline that is twice continuously differentiable at the 
secondary knots and continuous but with discontinuous derivatives at the primary 
knot locations t = 0, 1, 2, 3, 4. In terms of the parameterization (5.45), the knots and 
the smoothness at their locations are specified as the set 

{0, 0,0, 5,1, 1,1, 3,2, 2,2, 3,3, 3,3, 4,4,4, 4}. 


The knots are marked on the graph on the left in Fig. 5.122. 
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Fig. 5.122 Left: the knot locations on 0S. Right: the collocation points. 
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The 8 knots determine 16 functions b; ; for each component of W, generating a 
total of 32 basis functions. The graphs of the b;,; in parametric form are shown in 
Fig. 5.123. The first basis function has been modified to ensure periodic continuity 
at the knot locations t = 0, 4. 
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Fig. 5.123 The 16 B-spline basis functions b; ;[f], 0<t < 4. 
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To compute W, we must choose an appropriate number of collocation points, 
which can be placed anywhere because y is expected to be continuous on 0S. We 
select equally spaced points that avoid the primary knots. The set of 16 values of the 
parameter ¢ at our selection is 


ae eee 13 15 17 19 21 23 25 27 29 31 
8? 87? 8° 8° 8? 8° 8’? 87? 8? 8°? 8° 8° 8% 8’? 8° 8 
These values are marked on the graph on the right in Fig. 5.122. 

The same 16 collocation points are used for both components of W, yielding a 
system of 32 constraining equations that enables us to compute the coefficients Cg, j, j 
in (5.48). The 32 x 32 matrix of this system has a condition number of 245. 


5.12.4 Solution 


Approximate density. The two components of W are graphed in Figs. 5.124 and 
5.125. The latter also shows the spline knot locations. It can be seen that W agrees 
very well with the exact density w[x] = u[x], x € 0S, shown in Figs. 5.120 and 5.121. 


1 1 


Fig. 5.124 The components of [x] (Cartesian coordinates). 
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Fig. 5.125 The components of [x[t]] (parametric form). 
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Approximate solution. We use W and .% in (5.48) to compute an approximation 
ii to the exact solution u in S. The graphs of the two components of @ are shown in 


Fig. 5.126. 


Fig. 5.126 The components of [x], x € S. 


Error analysis. To visualize the difference between the approximate and exact 
densities, in Fig. 5.127 we graphed W — y relative to the maximum absolute value 
of yon ds. 


Fig. 5.127 The components of the relative error in W[x] (Cartesian coordinates). 


The same was done in parametric form in Fig. 5.128, which also shows the pri- 
mary and secondary knots and indicates that the relative error in the two components 
of W is about 0.05%. With such a small error on the boundary, we expect an equally 
good agreement in the interior of the domain. This is indeed confirmed by Fig. 
5.129, which displays the graphs of the components of #— u in S. 
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Fig. 5.128 The components of the relative error in W[x[t]] (parametric form). 


Fig. 5.129 The components of a[x] —ul[x], x € S. 


5.13 Remark. The set of coefficients {cj jj, c2;,;} in (5.48), computed with the 
B-spline basis shown in Fig. 5.123, is 


{{0.477, 0.491, 0.522, 0.631, 0.701, 1.84, 4.39, 8.07, 
10.3, 8.9, 6.62, 5.15, 4.53, 3.44, 1.83, 0.859}, 

{-14.8, -16.4, —20.8, —28.8, —34.2, —34.8, —36.7, —40.5, 
43., —36.4, 26.6, —21.1, —18.9, 17.7, —15.8, —15.}}. 
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5.13 Robin Problem in a Square Revisited: Piecewise Cubic 
Spline 


5.13.1 Summary 


This example illustrates the application of the alternative direct method to the Robin 
boundary value problem in a square domain. We compute the boundary condition 
from a known test solution, then use the collocation method with a piecewise cu- 
bic spline on equally spaced knots and collocation points to approximate the exact 
solution. It turns out that this alternative formulation is better conditioned than the 
version in Sect. 5.12. The results are validated through a comparison between the 
approximate and exact solutions. 


5.13.2 Problem Statement 


Domain boundary. The contour 0S is the square parameterized by 


t, 0<r<il, 0, O<r<l, 
1 1<t<2 -14+t, 1<t<2 
xlffJ=d 2” ee ee (5.49) 
3-t, 2<1<3, 1, 2<t<3, 
0, 3<t<4, 4-t 3<t<4. 


Its graph can be seen in Fig. 5.135. 

Governing equations. The displacement u satisfies the boundary value problem 
Zyoulx] =0, xES, 
(Tu) [x] + oufx] = # [x], x © OS. 


Test solution. The two components of the test solution u given by (5.3) are graphed 
in Fig. 5.130. 
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Fig. 5.130 The components of u[x], x € S. 
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Boundary data function. The function -% is given by 
(Tu) [x]+ oufx] = “ [x], xe OS, 


where, for computational purposes, we have chosen (without loss of generality) the 
particular 2 x 2 matrix 


The graphs of the two components of .% are displayed in Figs. 5.131 and 5.132, 
in Cartesian coordinates and parametric form, respectively. 
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Fig. 5.131 The components of .% [x] (Cartesian coordinates). 
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Fig. 5.132 The components of .% [x[t]] (parametric form). 
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5.13.3 Solution Procedure 


Method. We apply the alternative direct method, which reduces the problem to the 
boundary integral equation 


(Wo+ 51 G+V(oo) =V(%) ond, 
coded as 


f Plesloobldn+ f Dixso(eob)ar+} ols 
Topvpy TWeak([x] 


= $ D[x,y] op] db. (5.50) 
TWeak(x] 


This equation is solved numerically to obtain an approximate density @. In turn, @ 
is used in the representation 


ue] = f Dis.y]o(%—o@p)aty— J Pxslophldh, xeS, SSN) 
Tr Tr 


to generate an approximate solution @ in S. 


Density function. The exact density @ is the boundary stress vector Tu computed 
from the test solution u and given by (5.4); its components are graphed in Figs. 
5.133 and 5.134. 


i 1 


Fig. 5.133 The components of ¢[x] (Cartesian coordinates). 
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Fig. 5.134 The components of ¢[x{[f]] (parametric form). 


Numerical approximation. We compute @ by the collocation method with a 
B-spline basis of elements b;,;. Then the approximate density is sought in the form 


Lhe, (i,j) 
i ij 


@[x{¢]] = EEen, jbi le] } (5.52) 
ij 
where the coefficients Cg ;,; are determined by substituting (5.52) in (5.50). 

Since @ represents Tu on 0S and the latter is discontinuous at the corners, the bj, ; 
need to exhibit the same feature. In our case, we choose a piecewise cubic spline 
that is twice continuously differentiable at the secondary knots and discontinuous at 
the primary knot locations t = 0, 1, 2, 3, 4. In terms of the parameterization (5.49), 
the knots and the smoothness at their locations are specified as the set 


{0,0,0;5,1,.1¢1; 1,.2,.2, 2,2, 2,23, 3,3; 3;.454,4,4,4}. 


The knots are marked on the graph on the left in Fig. 5.135. 


Fig. 5.135 Left: the knot locations on 0S. Right: the collocation points. 
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The 8 knots determine 18 functions b; ; for each component of @, generating a 
total of 36 basis functions. Their graphs in parametric form are shown in Fig. 5.136. 
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Fig. 5.136 The 20 B-spline basis functions b;,;[t], 0<t < 4. 
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To compute @, we must choose an appropriate number of collocation points, 
which can be placed anywhere except at corners, because @ is expected to be dis- 
continuous there. Hence, we select collocation points that avoid the primary knots 
and are equally spaced. The set of 20 values of the parameter ¢ at our selection is 


1 3 17 9 I 13 3 17 19 21 23 5 27 29 31 33 7 37 39 
10° 10° 2? 10° 10° 10 10’ 2° 10’ 10’ 10’ 10° 27 10° 10° 10° 10° 2? 10’ 10S 


These values are marked on the graph on the right in Fig. 5.135. 

The same 20 collocation points are used for both components of @, yielding a 
system of 40 constraining equations that enable us to compute the coefficients cg ;j, ; 
in (5.52). The 40 x 36 matrix of this system has a condition number of 21. 


5.13.4 Solution 


Approximate density. The two components of @ are graphed in Figs. 5.137 and 
5.138. The latter also shows the locations of the primary and secondary knots. 


1 1 


Fig. 5.137 The components of [x] (Cartesian coordinates). 
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Fig. 5.138 The components of @[x[t]] (parametric form). 
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Approximate solution. We use @ and .% in (5.52) to compute an approximation 
ii to the exact solution u in S. The graphs of the two components of @ are shown in 


Fig. 5.139. 


Fig. 5.139 The components of di[x], x € S. 


Error analysis. To visualize the difference between the approximate and exact 
densities, in Fig. 5.140 we graphed @ — @ relative to the maximum absolute value 


of @ on OS. 


Fig. 5.140 The components of the relative error in @[x] (Cartesian coordinates). 


The same was done in parametric form in Fig. 5.141, which also shows the pri- 
mary and secondary knots and indicates that the relative error in the two components 
of @ is about 0.02%. With such a small error on the boundary, we expect an equally 
good agreement in the interior of the domain. 
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Fig. 5.141 The components of the relative error in @[x[t]] (parametric form). 


5.14 Remarks. (i) The amount of computation needed when we use (5.50) or (5.46) 
is about the same. This is expected because of the similarity in the structure of these 
equations. Also, the overall relative accuracy in both cases is similar (considering 
that the former calculates u and the latter calculates Tu on the boundary). The accu- 
racy for (5.50) is slightly better, which is to be expected because the B-spline basis 
for it has more elements. However, what is not expected is that the condition num- 
ber of the coefficient matrix resulting from the use of the (5.50) is 21, significantly 
smaller than the condition number of 245 when we use equation (5.46). Further 
work is required to see if this trend is confirmed and, thus, we are justified to use of 
the alternative direct method in preference to the ‘classical’ direct method. 

(ii) The set of coefficients {c; jj, c2;,;} in (5.52), computed with the B-spline 
basis shown in Fig. 5.136, is 


4£415,1;.18.6; 07.2, Ad, $77, 0.932), 3,19, 11.6,22.7, 99.4, 52, 
41.1, 25.3, —16.2, —12.5, —11.2, —8.45, —4.02, —0.707, 0.865}, 

£6.11, 7.33, 10., 13.9, 16.3, —52.2, —52., —52., —52.1, —52., —66.4, 
58.1, —45.2, —37., —33.7, 12.5, 13.2, 14.1, 14.7, 15.1} }. 


5.14 Dirichlet Problem in a Square: Classical Indirect Method 


5.14.1 Summary 


This example illustrates the application of the classical indirect method to the 
Dirichlet boundary value problem in a square domain. We compute the boundary 
condition from a known test solution, then use the collocation method with a piece- 
wise cubic spline on equally spaced knots and collocation points to approximate the 
exact solution. The results for the unknown density are validated by means of two 
different techniques. 
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5.14.2 Problem Statement 


Domain boundary. 0S is the square parameterized by 


t, 0<r<l, 0, 0<r<l, 
1, 1<t<2, t-1l, 1<t<2, 
x1[t] = — x2[t] = (5.53) 
3-t, 2<1<3, 1, Dots 3, 
0, 3<t<4, 4—-t 3<t<4. 


Its graph can be seen in Fig. 5.145. 


Governing equations. The displacement wu is the solution of the boundary value 
problem 


Z,oulx] =0, xES, 
ulx] = Aix], xe. 


Test solution. The two components of the test solution u given by (5.3) are graphed 
in Fig. 5.142. 


Fig. 5.142 The components of u[x], x € S. 


Boundary data function. The function Y is computed from the test solution u as 
Px} =ulx], xe OS. 


The graphs of its two components are displayed in Figs. 5.143 and 5.144, in Carte- 
sian coordinates and parametric form, respectively. 
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1 1 


Fig. 5.143 The components of #[x] (Cartesian coordinates). 


10 


Fig. 5.144 The components of Y[x{t]] (parametric form). 


5.14.3 Solution Procedure 


Method. We apply the classical indirect method, which reduces the problem to the 
boundary integral equation 


(Wo — 51) =P onds, 


coded as 


$ Plesloobldt— 4b] = Ab (5.54) 
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This equation is solved numerically to obtain an approximate density @. In turn, @ 
is used in the representation 


ulx] = $Phx,3] o@lyldy, xeSs (5.55) 
T 


to generate an approximate solution @ in S. 


Density function. The exact density @ does not represent either u or Tu on OS, 
which makes validation of the results somewhat more difficult. 


Numerical approximation. We compute @ by the collocation method with a 
B-spline basis of elements b; ;. Then the approximate density is sought in the form 


LL, i,j] 


@lx{t]] = TTe4jujlA (5.56) 
ij 


where the coefficients cq ;,; are determined by substituting (5.56) in (5.54). 

The construction of a spline basis requires that the smoothness of the spline at the 
knots must be the same as, or weaker than, that of the function being approximated. 
Therefore, since we do not know the behavior of @ at the corners, we have to work 
with a spline that is discontinuous at those four points. In our case, we choose a 
piecewise cubic spline that is twice continuously differentiable at the secondary 
knots and discontinuous at the primary knot locations t = 0, 1, 2, 3, 4. In terms of the 
parameterization (5.53), the knots and the smoothness at their locations are specified 
as the set 


{0, 0, 0,0, 5, 3,1, 1,1, 1, §, 3,2, 2,2, 2, 4, 2,3,3, 3,3, 2, 4,4,4,4, 4}. 


The knots are marked on the graph on the left in Fig. 5.145. 
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Fig. 5.145 Left: the knot locations on 0S. Right: the collocation points. 
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The 13 knots determine 24 functions b;,; for each component of @, generating 48 
basis functions. The graphs of the b;,; in parametric form are shown in Fig. 5.146. 


1 


1 


Fig. 5.146 The 48 B-spline basis functions b;,;[t], 0<t <4. 
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To compute @, we must choose an appropriate number of collocation points, 
which can be placed anywhere except at the corners t = 0, 1, 2, 3, 4, where @ might 
be discontinuous. Hence, we select collocation points that avoid the primary knots 
and are equally spaced. The set of 24 values of the parameter ¢ at our selection is 


eee: 73 11 13.5 17 19 7 23 
12° 4? 12? 12? 49 12? 12? 4? 12? 12? 4> 12? 


25 9 29 31 It 35 37 13 41 43 15 47 
12° 47 12° 12° 4° 12°? 12? 47 127 129 47 12 S° 


These values are marked on the graph on the right in Fig. 5.145. 

The same 24 points are used for both components of @, yielding a system of 48 
constraining equations that enable us to compute the numerical coefficients cg ;,; in 
(5.56). The 48 x 48 matrix of this system has a condition number of 41. 


5.14.4 Solution 


Approximate density. The two components of @ are graphed in Figs. 5.147 and 
5.148. These graphs strongly suggests that the density @ is continuous at the primary 
knot locations t = 0, 1, 2, 3, 4. 


1 1 


Fig. 5.147 The components of [x] (Cartesian coordinates). 
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Fig. 5.148 The components of @[x{t]] (parametric form). 
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Approximate solution. We use @ and # in (5.55) to compute an approximation a 
to the exact solution u in S. The graphs of the two components of #, shown in Fig. 
5.149, are virtually indistinguishable from those of u in Fig. 5.142. 


Fig. 5.149 The components of @i[x], x € S. 


Error analysis. A direct comparison between @ and ¢ is not possible because 
the latter is unknown. However, we can do this indirectly by using the fact that 
satisfies (see Sect. 5.14.3) 


(Wo-51I)e=f onods. 
The approximate density @ replaced in that equation produces an approximation Pp 
to FY. The difference Y — F relative to the maximum absolute value of Y on OS is 


graphed in Figs. 5.150 and 5.151. The latter shows that the relative error in the two 
components of # is less than 1%, which confirms the validity of our computation. 


0.01 


1 1 


Fig. 5.150 The components of the relative error in [x] (Cartesian coordinates). 
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Fig. 5.151 The components of the relative error in Y[x{¢]] (parametric form). 


A second way of comparing @ and @ is by looking at the difference between i 
in S, generated by (5.56) with @ replaced by @, and u in S, generated by the same 
formula with the exact density @. The difference #—u relative to the maximum 
absolute value of u in S, graphed in Fig. 5.152, shows good agreement between the 
computed and exact solutions, so, once again, we are confident that our numerical 
results are valid. 
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Fig. 5.152 The components of the relative error in [x], x € S. 


5.15 Remarks. (i) As mentioned earlier, the graphs of @ strongly suggest that the 
exact density @ is continuous at the corners of the domain, but that its derivatives 
are not. This feature should be taken into account when we construct a spline basis. 


(11) Computationally, the classical indirect method appears to be more sensitive 
to small errors than the direct method for similar domains with similar boundary 
conditions. As a result, the use of the former requires greater attention to be paid to 
computational accuracy. 
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(iii) The set of coefficients {c1;,;, c2,;,;} in (5.56), computed with the B-spline 
basis shown in Fig. 5.146, is 
{{—100., —102., —111., —133., —153., —164., —158., —96.8, 
— 20.4, 71., 162., 247., 259., 234., 198., 165., 
156., 154., 156., 131., 70., -22., -78.2, —101.}, 
{59 9A S, 197), — 9905 996. 307), A, 383 
133., —360., —374., —419., —455., —442., —347., —237., 
79.6, 58.4, —60.8, —52.2, —43.2, 42.9, —48.8, —54.1}}. 


5.15 Dirichlet Problem in an Asymmetric Domain: Classical 
Indirect Method 


5.15.1 Summary 


This example illustrates the application of the classical indirect method to the 
Dirichlet boundary value problem in a domain with a smooth asymmetric bound- 
ary. We compute the boundary condition from a known test solution, then use the 
collocation method with a piecewise linear spline on equally spaced knots and collo- 
cation points to approximate the exact solution. The results for the unknown density 
are validated indirectly by reference to the boundary data function. 


5.15.2 Problem Statement 


Domain boundary. The contour 0S is the curve parameterized by 


x1 [1] = 45 (7 + 8Cos[zt] + 4Cos[2z1]), 


<1<2. (5.57) 
x2[t] = sp (25 + 2Cos[3z1] + 20Sin[zz]), 


Its graph can be seen in Fig. 5.156. 


Governing equations. The displacement u is the solution of the boundary value 
problem 


Zyoulx] =0, xES, 


ulx]}= Pix], xeoSs. 


Test solution. The two components of the test solution u given by (5.3) are graphed 
in Fig. 5.153. 
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Fig. 5.153 The components of u[x], x € S. 
Boundary data function. The function Y is computed from the test solution u as 
Px] =ulx], xe OS. 


The graphs of its two components are displayed in Figs. 5.154 and 5.155, in Carte- 
sian coordinates and parametric form, respectively. 


Fig. 5.154 The components of # [x] (Cartesian coordinates). 
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Fig. 5.155 The components of #[x[t]] (parametric form). 
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5.15.3 Solution Procedure 


Method. We apply the classical indirect method, which reduces the problem to the 
boundary integral equation 
(Wo-51Ie=Ff onds, 


coded as 


¢ Plyo@bldty—4ob4 = A (5.58) 
Topvix| 


This equation is solved numerically to obtain an approximate density @. In turn, @ 
is used in the representation 


ulx] = f Plxylo@bldh;, xes (5.59) 
T 


to generate an approximate solution @ in S. 


Density function. The exact density @ does not represent either u or Tu on 0S, 
which makes validation of the results somewhat more difficult. 


Numerical approximation. We compute @ by the collocation method with a 
B-spline basis of elements b;,;. Then the approximate density is sought in the form 


rXe1 4,51, 
ij 


@[x{¢]] = ee nth (5.60) 
ij 
where the coefficients cy ;,; are determined by substituting (5.60) in (5.58). 

The construction of a spline basis requires the smoothness of the spline at the 
knots to be the same as, or weaker than, that of the approximated function. Since 
the parameterization of 0S is smooth, we expect the exact density @ to be smooth on 
OS as well, so we construct a spline with the same feature. In our case, we choose a 
piecewise linear spline continuous at all the knots. In terms of the parameterization 
(5.57), the knots and the smoothness at their locations are specified as the set 


{0, toa dp pipe pep ip ail 


15 8 17 9 19 10 3 I 23 12 25 13, 27,9} 
147.7) 14> 7) 14> 7929 7? 149 7? 149 7? 149 “S° 


The knots are marked on the graph on the left in Fig. 5.156. 
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Fig. 5.156 Left: the knot locations on 0S. Right: the collocation points. 
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The 28 knots determine 28 functions 5; ; for each component of @, generating 56 
basis functions whose graphs, in parametric form, are shown in Fig. 5.157. 
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Fig. 5.157 The 28 B-spline basis functions b;,;[t], 0 <t < 2. 
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To compute @, we must choose an appropriate number of collocation points. The 
expected smoothness of @ allows us to place these points almost anywhere. The set 
of 28 equally-spaced values of the parameter ¢ at our selection is 
2149 792? 7? 14? 7? 14° 7? 14? 


19 10 3 Il 23 12 25 13 27 
0749792979 149 7? 149 7? 4S 


g90t.1.3 2.5 3 14 9 5 H 6 By 
{0, >a G4 ’ 
17 9 
14° 7 


These values are marked on the graph on the right in Fig. 5.156. 

The same 28 points are used for both components of @, yielding a system of 56 
constraining equations that enable us to compute the coefficients cg; ; in (5.60). The 
56 x 56 matrix of this system has a condition number of 27. 


5.15.4 Solution 


Approximate density. The components of @ are graphed in Figs. 5.158 and 5.159. 
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Fig. 5.158 The components of [x] (Cartesian coordinates). 
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Fig. 5.159 The components of @[x[t]] (parametric form). 
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Error analysis. A direct comparison between @ and @ is not possible because the 
latter is unknown. However, we can do this indirectly by using the fact that @ is the 
exact solution of the equation 


(Wo-—51)e=f onds. 


The approximate density @ replaced in that equation produces an approximation Y 
to Y. Both functions Y and Y are graphed in Figs. 5.160 and 5.161. 


Fig. 5.160 The components of “[x] (orange line) and Ax] (blue line) in Cartesian coordinates. 


10 
0! 4 
0 1 2 
0} 1 
1 2 
—80 


Fig. 5.161 The components of A [x[r]] (orange line) and Y|[x[¢]] (blue line) in parametric form. 


We note that the agreement between the two is not very good. The first compo- 
nent of Y in Fig. 5.160 is twisted on the right-hand side of the domain, which is 
also where the parameterization begins and ends (see Fig. 5.161). Such behavior is 
characteristic of the errors produced by the classical indirect method and can be- 
come quite severe. We expected this to happen because we chose equal spacing of 
the parameter values for both the spline knots and collocation points, and not for 
their geometric positions to the boundary. As Fig. 5.156 shows, the density of these 
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points is high on the left segment of 0S and low on the right one. Subsequent calcu- 
lations, not shown here, show that a more uniform redistribution of the spline knots 
and collocation points along the boundary rather than with respect to their parameter 
values yields significantly better results. 


5.16 Remarks. (i) In the classical indirect method, the collocation points must 
be chosen carefully. Numerical experimentation shows that placing them midway 
between the spline knots when using piecewise linear splines produces an ill- 
conditioned problem. To avoid this, in our example we chose the collocation points 
at the knot locations. 


(11) In Fig. 5.157, the first basis function was modified to exhibit periodicity and 
thus be consistent with the expected continuity of @. 


(iii) The set of coefficients {c1;,;, c2,;,;} in (5.60), computed with the B-spline 
basis shown in Fig. 5.157, is 


{{167., 359., 416., 358., 267., 203., 175., 160., 150., 148., 147., 
135., 108., 68., 23.6, —13.9, —40., —56.6, —67.7, —76.6, —88., 
105., —126., —157., —203., —242., —214., —65.3}, {—1760., 
1580., —1150., —709., —391., —204., —94.8, —37.8, —24.3, 
25.4, 22.5, -17.5, —15.4, -15.4, 14.1, —8.68, 1.82, 14., 22.6, 
24., 14., -19.4, 89.1, —201., -387., —699., —1130., —1570.}}. 


5.16 Neumann Problem in an Ellipse: Classical Indirect Method 


5.16.1 Summary 


This example illustrates the application of the classical indirect method to a Neu- 
mann problem in a domain with a smooth boundary. We compute the boundary 
condition from a known test solution, then use the collocation method with a piece- 
wise cubic spline on equally spaced knots and collocation points to approximate 
the exact density. We also discuss the occurrence of ill-conditioning and indicate a 
technique for mitigating its effect. Since the exact density is unknown, the results 
are validated through an indirect procedure. 


5.16.2 Problem Statement 


Domain boundary. 0S is the ellipse parameterized by 


x1[t]=14Cos[at], x2[f] = 5Sin[ae]+5, O<1<2. (5.61) 


Its graph can be seen in Fig. 5.165. 
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Governing equations. The displacement u satisfies the boundary value problem 
Zyoulx]=0, xES, 
(Tu) [x] = T,oulx4] = Qx], xe OS. 


Test solution. The two components of the test solution u given by (5.3) are graphed 
in Fig. 5.162. 
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Fig. 5.162 The components of u[x], x € S. 


Boundary data function. The function 2 = Tu is given by (5.4). The graphs of 


its two components are displayed in Figs. 5.163 and 5.164, in Cartesian coordinates 
and parametric form, respectively. 
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Fig. 5.164 The components of 2[x[f]] (parametric form). 
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5.16.3 Solution Procedure 


Method. We apply the classical indirect method, which reduces the problem to the 
boundary integral equation 


(Wy +31) w=2 onds, 
coded as 


¢ (KeDb.y)owblah +} vb = 2 (5.62) 
Topvpy 
This equation is solved numerically to obtain an approximate density wy. In turn, w 
is used in the representation 


ule] = f Dixy owly]dh, xeS (5.63) 
Tr 


to generate an approximate solution @ in S. 


Density function. The exact density y does not represent either u or Tu on OS, 
which makes validation of the results somewhat more difficult. 


Numerical approximation. We compute by the collocation method with a 
B-spline basis of elements b;,;. Then the approximate density is sought in the form 


XXe14,j51,j[¢] 


Wx([e\] = Sea bisll ; (5.64) 
ij 


where the coefficients cy ;,; are determined by substituting (5.64) in (5.62). 

The construction of a spline basis requires the smoothness of the spline at the 
knots to be the same as, or weaker than, that of the approximated function. Since the 
parameterization (5.61) is smooth, we expect y to be equally smooth on 0S, so we 
may construct an approximate density W with any degree of smoothness. Therefore, 
to ensure maximum smoothness, we choose a piecewise cubic spline twice continu- 
ously differentiable at the knots and modified to enforce periodicity and continuity 
at the knots corresponding to t = 0, 2. In terms of the parameterization (5.61), the 
knots and the smoothness at their locations are specified as the set 


11125 7 43 5 I 
{0 53 5559 39 89 ds 39 39 Fs Gr ZF 


The knots are marked on the graph on the left in Fig. 5.165. 
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Fig. 5.165 Left: the knot locations on 0S. Right: the collocation points. 


202 5 Computational Examples 


The 12 knots determine 12 functions b; ; for each component of , generating a 
total of 24 basis functions. The graphs of the b;,; in parametric form are shown in 
Fig. 5.166. 
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Fig. 5.166 The 12 B-spline basis functions b;, ;[t], 0<t < 2. 


To compute Y, we must choose an appropriate number of collocation points. 
Owing to the smoothness of 0S, these points can be placed almost anywhere. How- 
ever, computational experiments (not shown here) indicate that the typical choice 
for these points generates a singular coefficient matrix. Consequently, to help with 
the computational stability of the method, we have augmented the original set with 
additional locations, which gives rise to an overdetermined system. We select col- 
location points that avoid the knots and are equally spaced. The set of 36 values of 
the parameter f at our selection is 

{A 1151735 UB7573 197 UL 23 
24? 12? 8? 24? 4? 24? 8? 12? 24 24? 12 8? 24? 4? 24? 8? 12? 24? 


25 13 9 29 5 31 Il 17 35 37 19 13 41 7 43 15 23 47 
24> 12? 8) 24 4) 24> 8 12> 24> 24? 12> B? 24? 49 24> BS 129 WS- 


These values are marked on the graph on the right in Fig. 5.165. 
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The same 36 points are used for both components of the approximation W, yield- 
ing a system of 72 constraining equations that enable us to compute the coefficients 
Co,i,j im (5.64). The 72 x 24 matrix of this overdetermined system has a condition 
number of 12. 


5.16.4 Solution 


Approximate density. The two components of W are graphed in Figs. 5.167 and 
5.168. These graphs strongly suggests that our conjecture regarding the continuity 
of the exact density y on 0S is correct. 


Fig. 5.167 The components of [x] (Cartesian coordinates). 
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—1000 
Fig. 5.168 The components of W[x[t]] (parametric form). 
Approximate solution. We use W in (5.63) to compute an approximation i to the 


exact solution u in S. The graphs of the components of # are shown in Fig. 5.169. 
There is virtually no visual distinction between # in Fig. 5.169 and u in Fig. 5.162. 


204 5 Computational Examples 


—100 


0 


2 


N 


Fig. 5.169 The components of [x], x € S. 


Error analysis. A direct comparison between W and y is not possible because the 
latter is unknown. However, we can do this indirectly by using the fact that u is 
computed from y by means of (5.63). The approximate density W replaced in the 
restriction of this formula to the boundary produces an approximation i to u on OS. 
The difference ii[x] — u[x], x € OS, relative to the maximum absolute value of u on 
OS is graphed in Figs. 5.170 and 5.171. The latter shows that the relative error in the 
two components of # on the boundary is less than 1%, which confirms the validity 
of our computation. 
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Fig. 5.171 The components of the relative error in fi[x{t]], 0 <t <2 (parametric form). 
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5.17 Remarks. (i) Computationally, the classical indirect method appears to be 
more sensitive to small errors than the direct method for similar domains with simi- 
lar boundary conditions. As a result, the use of the former requires greater attention 
to be paid to computational accuracy. Here, we used surplus collocation points to 
diminish the instability inherent in this technique. 

(ii) The set of coefficients {c1 jj, c2;,;} in (5.64), computed with the B-spline 
basis shown in Fig. 5.166, is 


{{133., 60., 42.7, 20., 19.2, 10.6, 9.37, 1.19, 11.7, 10., 
39.2, 61.6}, { -1200., —734., -377., —207., —116., 
81.2, —63., 37.7, —53.4, —88.5, —220., —494.}}. 


5.17 Neumann Problem in a Square: Classical Indirect Method 


5.17.1 Summary 


This example illustrates the application of the classical indirect method to a Neu- 
mann problem in a square domain, where the boundary condition is computed from 
a known test solution. We use the collocation method with a piecewise cubic spline 
on unequally spaced knots and collocation points to construct an approximation of 
the original solution from the computed boundary data function. The validation pro- 
cess, carried out through a comparison between the exact and approximate solutions, 
indicates that the classical indirect method is not recommended for use in domains 
with corners. 


5.17.2 Problem Statement 


Domain boundary. 0S is the square parameterized by 


t, 0<r<il, 0, 0<r<il, 
1 l<t<2 t-1, 1<t<2 
xI[fJ= eo” Set = =! (5.65) 
3-t, 2<t<3, 1, 2<t<3, 
0, 3<t<4, 4—-t 3<t<4. 


Its graph can be seen in Fig. 5.175. 
Governing equations. The displacement wu is the solution of the boundary value 
problem 

Zyoulx] =0, xES, 

(Tu) [x] = T.oulx] = Qix], xe OS. 


5 Computational Examples 


206 
Test solution. The two components of u given by (5.3) are graphed in Fig. 5.172. 


Fig. 5.172 The components of u[x], x € S. 


Boundary data function. The function 2 = Tu is given by (5.4). The graphs of 
its two components are displayed in Figs. 5.173 and 5.174, in Cartesian coordinates 


and parametric form, respectively. 


1 


1 
Fig. 5.173 The components of 2[x] (Cartesian coordinates). 
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Fig. 5.174 The components of 2[x[t]] (parametric form). 
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5.17.3 Solution Procedure 


Method. We apply the classical indirect method, which reduces the problem to the 
boundary integral equation 


(Wo +5)w=2 onds, 
coded as 


¢ (LoDkxy) owblahy +} ved = 2b (5.66) 


TCPV ix] 


This equation is solved numerically to obtain an approximate density w. In turn, w 
is used in the representation 


u[x] = $ Dis.y owl|dI,, xEs (5.67) 
F 


to generate an approximate solution d in S. 


Density function. The exact density y does not represent either u or Tu on OS, 
which makes validation of the results somewhat more difficult. 


Numerical approximation. We compute by the collocation method with a 
B-spline basis of elements b; ;. Consequently, the approximate density is sought 
in the form 


LL. ji, le) 
fs ig 


Wlx[e]] = LDe24, jij ’ (5.68) 
ij 


where the numerical coefficients cq;,; are determined by substituting (5.68) in 
(5.66). 

The construction of a spline basis requires that the smoothness of the spline at the 
knots must be the same as, or weaker than, that of the function being approximated. 
Therefore, since we do not know the behavior of y at the corners, we have to work 
with a spline that is discontinuous at those four points. In this case, we choose a 
piecewise cubic spline that is twice continuously differentiable at the secondary 
knots and discontinuous at the primary knot locations ¢ = 0, 1, 2, 3, 4. 

According to the calculations below, y appears to have a singularity at the cor- 
ners, so the collocation method would perform very poorly if we were to use a typ- 
ical uniform distribution of spline knots. Instead, we choose knots and collocation 
points that are densely clustered about the parametric values t = 0, 1, 2, 3, 4 corre- 
sponding to the four corners of the boundary. This type of choice helps us judge the 
behavior of y at those points. In terms of the parameterization (5.65), the knots and 
the smoothness at their locations are specified as the set 
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{0, 0, 0, 0, 0.00302, 0.0141, 0.107, 0.5, 0.893, 0.986, 0.997, 1, 1, 1, 1, 1.003, 
1.014, 1.107, 1.5, 1.893, 1.986, 1.997, 2, 2, 2, 2, 2.003, 2.014, 2.107, 2.5, 
2.893, 2.986, 2.997, 3, 3, 3, 3, 3.003, 3.014, 3.107, 3.5, 3.893, 3.986, 3.997, 
4,4, 4, 4}. 


The knots are marked on the graph on the left in Fig. 5.175. We draw attention to the 
fact that, owing to the closeness of the knots near the corners, 16 of the 28 secondary 
knot locations are hidden under the large points designating the four primary knot 
locations. 


—— 
T 
oo 


Fig. 5.175 Left: the knot locations on 0S. Right: the collocation points. 


The 32 knots determine 44 functions b;,; for each component of W, generating a 
total of 88 basis functions. The graphs of the b;,; in parametric form are shown in 
Fig. 5.176. Since the visual aspect of the functions constructed in the vicinity of the 
corners is not very helpful, in Fig. 5.177 we have redrawn the restrictions of the first 
four to the interval 0 < t < 0.017, indicating the three spline knot locations and the 
four collocation points on them by black and white circular markers, respectively. 

To compute W, we must choose an appropriate number of collocation points. The 
set of 64 unequally spaced values of the parameter ¢ at our selection is 


{0.001007, 0.002014, 0.006725, 0.01043, 0.04506, 0.07599, 0.2379, 
0.3690, 0.6310, 0.7621, 0.9240, 0.9549, 0.9896, 0.9933, 0.9980, 
0.9990, 1.001, 1.002, 1.007, 1.010, 1.045, 1.076, 1.238, 1.369, 1.631, 
1.762, 1.924, 1.955, 1.990, 1.993, 1.998, 1.999, 2.001, 2.002, 2.007, 
2.010, 2.045, 2.076, 2.238, 2.369, 2.631, 2.762, 2.924, 2.955, 2.990, 
2.993, 2.998, 2.999, 3.001, 3.002, 3.007, 3.010, 3.045, 3.076, 3.238, 
3.369, 3.631, 3.762, 3.924, 3.955, 3.990, 3.993, 3.998, 3.999}. 


5.17 Neumann Problem in a Square: Classical Indirect Method 
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Fig. 5.176 The 44 B-spline basis functions b; ;[f], 0<t < 4. 
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Fig. 5.177 The first four B-spline basis functions b;, ;[t], 0 < t < 0.017. 


These values are marked on the graph on the right in Fig. 5.175. As in the case of 
the knots, only 32 of the 64 collocation points are visible, the remaining 32 being 
hidden under the four white corner markers. However, a few of the latter can be seen 
in Fig. 5.177; they correspond to the parametric locations 


t = 0.001007, 0.002014, 0.006725, 0.01043. 


The same 64 collocation points are used for both components of W, yielding a 
system of 128 constraining equations that enables us to compute the coefficients 
Co,i,j in (5.68). The 88 x 128 matrix of this overdetermined system has a condition 
number of 281. 


5.17.4 Solution 


Approximate density. The two components of W/[x] are shown in Figs. 5.178 and 
5.179. These graphs strongly suggest that the density y[x/t]] is discontinuous at the 
primary knot locations tf = 0, 1, 2,3, 4. Additional computation shows that when 
supplementary spline knots and collocation points are placed closer to the corners, 
the spike in W[x[t]] at these locations increases in height, which supports the conjec- 
ture that y is singular at the four corners. 


1 


Fig. 5.178 The components of [x] (Cartesian coordinates). 
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Fig. 5.179 The components of [x[t]] (parametric form). 


Error analysis. A direct comparison between W and y is not possible because the 
latter is unknown. However, we can do this indirectly by using the fact that u is 
computed from y by means of (5.67). The approximate density W replaced in the 
restriction of this formula to the boundary produces an approximation i to u on OS. 
The graphs of the components of both uw (blue line) and # (magenta line) for x € 0S 
are shown in Fig. 5.180. 


1 


Fig. 5.180 The components of u[x] (blue line) and a[x] (magenta line), x € OS. 


These graphs tell us that the approximate and exact densities are in at least partial 
agreement for this example. However, the approximate density is not very accurate. 
Greater accuracy can, of course, be achieved by means of additional spline knots and 
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collocation points that are even more densely packed around the corner locations. 
But this procedure becomes increasingly ill behaved. Here, we have used t values 
within 0.001 of a corner. A computation attempting to reduce this distance to 0.0001 
requires a floating-point accuracy of over 70 digits, which very rapidly makes the 
problem unsolvable. Consequently, we do not recommend this method for problems 
involving boundaries with corners. 


5.18 Remarks. (i) As illustrated in Sect. 5.16, the classical indirect method is com- 
putationally acceptable for smooth boundary curves. 


(ii) The set of coefficients {c;,;, c2;,;} in (5.68), computed with the B-spline 
basis shown in Fig. 5.176, is 


{{14., 16.2, 16.2, 19.9, 22.4, 54.6, —40.6, —10., 9.36, 33.7, 78.1, 178., 
3.09, 25.9, 12.1, 5.73, 2.75, 11.1, 17.1, —3.28, 2.4, 7.27, 16.7, 38.7}, 
{-86.8, —57.4, —56.7, 77.1, —135., —268., —205., —169., —159., 
175., —219., —286., —470., —264., —167., —124., -117., 
156., —79.3, —51.2, —35., 29.1, —33.4, —42.6}}. 


5.18 Robin Problem in An Ellipse: Classical Indirect Method 


5.18.1 Summary 


This example illustrates the application of the classical indirect method to a Neu- 
mann problem in a domain with a smooth boundary. We compute the boundary 
condition from a known test solution, then use the collocation method with a piece- 
wise cubic spline on equally spaced knots and collocation points to approximate the 
exact density. Since the latter is unknown, the results are validated through an indi- 
rect procedure. We also make a comparison between the classical indirect method 
for this problem and for the Neumann problem discussed in Sect. 5.16. 


5.18.2 Problem Statement 


Domain boundary. QS is the ellipse parameterized by 
x1[f] =1+Cos[at], x2[t] = 4Sin[at]+ 3, O<t<2. (5.69) 


Its graph can be seen in Fig. 5.184. 


Governing equations. The displacement wu is the solution of the boundary value 
problem 
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Z,oulx] =0, xES, 
(Tu) [x] + oulx] = 4 |x], x€ OS. 


Test solution. The two components of the test solution u given by (5.3) are graphed 
in Fig. 5.181. 


Fig. 5.181 The components of u[x], x € S. 


Boundary data function. The function -% is given by 
(Tu) [x] + oulx] = 4 |x], x€ OS, 


where, for computational purposes, we have chosen the particular 2 x 2 matrix 


The graphs of the two components of .% are displayed in Figs. 5.182 and 5.183, 
in Cartesian coordinates and parametric form, respectively. 
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Fig. 5.182 The components of .% [x] (Cartesian coordinates). 
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Fig. 5.183 The components of .% [x[t]] (parametric form). 


5.18.3 Solution Procedure 


Method. We apply the classical indirect method, which reduces the problem to the 
boundary integral equation 


(Ws +3)9+Vo(o9) =H onds, 
coded as 
¢ (LoDkx.»)ooblah+hel+ ¢ Dkx.ylo(ogblar = bx, 


TCpv ix] TWeak{x] 
(5.70) 


This equation is solved numerically to obtain an approximate density @. In turn, @ 
is used in the representation 


u[x] = f Dis.y oglyldIy, xES (5.71) 
r 


to generate an approximate solution @ in S. 


Density function. The density @ does not represent either u or Tu on 0S, which 
makes validation of the results somewhat more difficult. 


Numerical approximation. We compute @ by the collocation method with a 
B-spline basis of elements 5; ;. Then the approximate density is sought in the form 


Le i,;51,;le] 
~ ij 


§[x{¢]] = Ld eo, ibis] (5.72) 
ij 


where the coefficients cq ;,; are determined by substituting (5.72) in (5.70). 
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The construction of a spline basis requires that the smoothness of the spline at the 
knots must be the same as, or lower than, that of the function being approximated. 
Since the parameterization (5.69) is smooth, we expect @ to be equally smooth on 
OS, so we may construct an approximate density @ with any degree of smoothness. 
Therefore, to ensure maximum smoothness, we choose a piecewise cubic spline that 
is twice continuously differentiable at the knots and has been modified to enforce 
periodicity and continuity at the knots corresponding to tf = 0, 2. In terms of the 
parameterization (5.69), the knots and the smoothness at their locations are specified 
as the set 


Fig. 5.184 Left: the knot locations on 0S. Right: the collocation points. 


The 12 knots determine 12 functions b; ; for each component of , generating a 
total of 24 basis functions. The graphs of the b;,; in parametric form are shown in 
Fig. 5.185. The first, second, and last basis functions have been modified to ensure 
periodic continuity at tf = 0, 2. 

To compute @, we must choose an appropriate number of collocation points. 
Owing to the smoothness of 0S, these points can be placed almost anywhere. How- 
ever, computational experiments (not shown here) indicate that the typical choice for 
these points generates a singular coefficient matrix. Consequently, to help with the 
computational stability of the method, we have augmented our original set of cen- 
trally placed locations with additional ones, which gives rise to an overdetermined 
system. The set of 36 equally spaced values of the parameter f at our selection is 

(A 1151735 UB7573 197 UL 23 
24? 12? 8? 24? 4? 24? 8? 12? 24 24 12? 8? 24? 4? 24? 8? 12? 24? 


25 13 9 29 5 31 Il 17 35 37 19 13 41 7 43 15 23 47 
24> 12> 8> 24> 49 24> B? 12) 24> 24? 12? B? 24? 4 24> B? 12? 4S 


These values are marked on the graph on the right in Fig. 5.184. 

The same 36 collocation points are used for both components of W, yielding a 
system of 72 constraining equations that enables us to compute the coefficients Cg, ;j, j 
in (5.72). The 72 x 24 matrix of the overdetermined system has a condition number 
of 8. 
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Fig. 5.185 The 12 B-spline basis functions b;,;[t], 0 <t < 2. 


5.18.4 Solution 


Approximate density. The components of @ are graphed in Figs. 5.186 and 5.187. 


2 


Fig. 5.186 The components of [x] (Cartesian coordinates). 
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Fig. 5.187 The components of @[x/t]] (parametric form). 


Approximate solution. We use @ in (5.71) to compute an approximation @ to the 
exact solution u in S. The graphs of the two components of # are shown in Fig. 
5.188. There is virtually no visual distinction between i# in these graphs and u in 
Fig. 5.181. 


2 2 


Fig. 5.188 The components of @[x], x € S. 


Error analysis. A direct comparison between @ and @ is not possible because the 
latter is unknown. However, we can do this indirectly by using the fact that u is 
computed from @ by means of (5.71). The approximate density W replaced in the 
restriction of this formula to the boundary produces an approximation i to u on OS. 
The difference [x] — u[x], x € 0S, relative to the maximum absolute value of u on 
OS is graphed in Figs. 5.189 and 5.190. The latter shows that the relative error in the 
two components of i on the boundary is less than 1%, which confirms the validity 
of our computation. 

Another procedure for illustrating the validity of our method is to compare the 
result for the Robin problem obtained above with that for the Neumann problem in 
Sect. 5.16 since in both cases the test solution u in S has the same representation; 
that is, equations (5.63) and (5.71) are essentially the same and generate the same 
function u, which means that the exact densities @ and y must coincide. 


218 5 Computational Examples 


1/2 1 3/2 3 
0.01 } 


Fig. 5.190 The components of the relative error in @[x[t]], 0 << <2 (parametric form). 


However, their approximations @ and W do not coincide because their compu- 
tational constructions are different, even though they are based on boundary data 
functions 2 and .% in (5.64) and (5.72) that are both computed from the test solu- 
tion (5.3). The difference @ — y relative to the maximum absolute value of y on 0S 
is graphed in Fig. 5.191. These graphs show that the relative size of the difference is 
less than 0.5%, which validates the results of the classical indirect method for both 
the Neumann and Robin boundary value problems. 
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Fig. 5.191 The components of the relative difference @[x[t]] — w[x[r]]. 
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5.19 Remarks. (i) We mention that the density @ in this problem suffers from the 
same computational issues as the density y in Sect. 5.17 as they are both used in the 
same integral representation of the exact solution u. Consequently, if the boundary 
has corners, we would encounter the same computational difficulties when con- 
structing @ as we did for w. 

(ii) The set of coefficients {c1 jj, c2;,;} in (5.72), computed with the B-spline 
basis shown in Fig. 5.185, is 


{{133., 61.4, 42.8, 20.4, 19.5, 10.5, 9.74, 0.881, 11.4, 


10.3, 37.9, 61.5}, {—1200., —735., —377., —208., —116., 
81.1, —62.7, —37.1, —53.9, —88., —220., —496.}}. 


5.19 Mixed Boundary Conditions: Piecewise Cubic Spline 


5.19.1 Summary 


This example illustrates the application of the direct method to a mixed problem 
in a semicircular domain with Dirichlet conditions prescribed on the diameter and 
Neumann conditions on the arc. Both data functions are computed from a known test 
solution. We use the collocation method with a piecewise cubic spline on equally 
spaced knots and collocation points to construct an approximation of the original 
solution from the exact boundary data functions. The results are validated through a 
comparison between the exact and approximate solutions. 


5.19.2 Problem Statement 


Domain boundary. 0S is a closed curve consisting of a straight line segment 0S, 
and the upper half circle 0S that has this segment as its diameter. It is parameterized 
by 


2r, 0<t<1, 
x1[¢] = 
1+Cos[z(t—-1)], 1<t<2, 
(5.73) 
<t< 
x2] = . 0<t<1, 
Sin[x(¢—1)], 1<1t<2, 


and its graph can be seen in Fig. 5.201. 
Governing equations. The displacement u is the solution of the boundary value 
problem 

Zyoulx] =0, xES, 

ulx] = Pix], xe dS), 

(Tu) [x] = Qa], x EAS). 
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Test solution. The two components of the test solution u given by (5.3) are graphed 
in Fig. 5.192. 
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Fig. 5.192 The components of u[x], x € S. 


Boundary data functions. The function Y is given by (5.3) on 0S,; its compo- 
nents are graphed in Fig. 5.193 and 5.194, in Cartesian coordinates and parametric 
form, respectively. 


Fig. 5.193 The components of Y |x], x € S$, (Cartesian coordinates). 
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Fig. 5.194 The components of Y[x{t]], 0 <t < 1 (parametric form). 
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The function 2 is given by (5.4) on AS»; its components are displayed in Figs. 
5.195 and 5.196. 


Fig. 5.195 The components of 2[x], x € 0S2 (Cartesian coordinates). 
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Fig. 5.196 The components of 2[x[t]], 1 <t <2 (parametric form). 


5.19.3 Solution Procedure 


Method. We apply the direct method, which reduces the problem to the boundary 
integral equations 


Vo(@) = (Wo + sl)P on 0S}, 


codes as 


¢ — Dixsoobldh=f — PyloPbldh+4 Abi, xe a5, 
TWeak{x] TCPV ix] 
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where @[x] = (Tu) [x], and 
(Wo+51)Q =Vo>2 ondSp, 


coded as 


¢  Plxvlo@bldh+sekl=p — DixyJoQpldh, xe a, 
Topvpy TWeak(x] 
(5.75) 


where @[x] = u[x], x € AS2. These equations are solved numerically to obtain an 
approximate density @. In turn, @ is used in the representation formula 


ud=f Drlooblah+ f Dixsle2blah 


-$ PkyloPbldy—¢ Pleslo@bldh, xes, (5.76) 
Os, : Sy s 


to generate an approximate solution d@ in the domain S. 


Density function. The exact density @ on 0S} is the restriction to 0S» of the test 
solution u; its components are graphed in Figs. 5.197 and 5.198 in Cartesian coor- 
dinates and parametric form, respectively. 
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Fig. 5.197 The components of g[x], x € 0S2 (Cartesian coordinates). 
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Fig. 5.198 The components of @[x[t]], 1 <+ <2 (parametric form). 
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The exact density @ on 0S is the restriction to 0S, of the boundary stress vector 
Tu computed from the test solution and given by (5.4); its components are graphed 
in Figs. 5.199 and 5.200. 
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Fig. 5.199 The components of @[x], x € 0S; (Cartesian coordinates). 
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Fig. 5.200 The components of @[x(t]], 0 <+¢< 1 (parametric form). 


Numerical approximation. We compute @ by the collocation method with a 
B-spline basis of elements b; ;. Then the approximate density is sought in the form 


Vdc jbi,jlt] 
erty _ | td 
Olx{[t]] = Ldc2 5, jbi, jt] 
j 


i 


(5.77) 


where the numerical coefficients cq ;,; are determined by substituting (5.77) in 
(5.74) and (5.75). 

Since @ is discontinuous at the two corners of 0S because of the boundary con- 
dition switch between Y and J at those points, the b;,; need to exhibit the same 
feature. In this case, we choose a piecewise cubic spline that is twice continuously 
differentiable at the secondary knots and discontinuous at the primary knot locations 
t = 0, 1, 2. In terms of the parameterization (5.73), the knots and the smoothness at 
their locations are specified as the set 


1 ti2.5 743 5 U 
{0,0,0,0,%,9,35.42 11,11, 2 3, 3,3, Z, 2,2, 2,2}. 
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The knots are marked on the graph on the left in Fig. 5.201. 
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Fig. 5.201 Left: the knot locations on 0S. Right: the collocation points. 


The 12 knots determine 18 functions b;,; for each component of W, generating a 
total of 36 basis functions. The first nine have support only on 0S, and the remain- 
ing nine have support only on OS. The graphs of the b; ; in parametric form are 
shown in Fig. 5.202. 

The computed density @[x|[t]] that approximates (Tu)|[x[t]] on 0 <1 < 1 is ex- 
pected to ‘join’ the boundary condition function 2[x{t]] (defined on 1 <t < 2) 
with jump discontinuities at t = 0, 1, 2 because the parameterized boundary con- 
tour 0S = 0S; US» has corners at these locations. Therefore, the first nine basis 
functions b; ; for 0S, in Fig. 5.202 are acceptable in their current form. On the other 
hand, the computed density @[x|t]] that approximates u[x/t]] on 1 <t <2 is expected 
to join the boundary condition function Y|x[t]] (defined on 0 < t < 1) without a 
jump at those locations, which suggests that we should impose an additional re- 
striction on the second nine basis functions to guarantee continuity at tr = 0, 1, 2. 
Specifically, we choose the coefficients at positions 9 and 18 in Fig. 5.202 to be 


@[x{1]] = u[x[1]] = {1.85013, 95.3664}, 
= u|x(0]] = {0.476909, — 14.8369}. 


This means that the boundary element method will compute only seven of the nine 
basis functions for each of the two components on 0S). 

To compute @, we must choose an appropriate number of collocation points, 
which can be placed anywhere on the boundary except at the corners, where @ is 
not defined. The set of values of the parameter f¢ at our selection is 


{1 15 7 1 1 13.5 17 19 7 23 25 3 29 31 I 35 
18> 6? 18> 18? 2° 18> 18> 6 18? 18> 6 18> 18> 2» 18> 18> 6» 8S" 


These values are marked on the graph on the right in Fig. 5.201. 

The first nine collocation points lie on 0S; and are used to determine the coeffi- 
cients of the first nine basis functions for x € 0S,; the other nine lie on 0S and are 
used to determine the coefficients of the last seven basis functions for x € Sz. The 
same 18 collocation points are used for both components of @, yielding a system 
of 36 constraining equations that enables us to compute the numerical coefficients 
Cq,i,j i (5.77). The overdetermined 36 x 32 matrix of this system has a condition 
number of 428. 
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Fig. 5.202 The 18 B-spline basis functions b;, ;[t], 0 <t < 2. 
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5.19.4 Solution 


Approximate density. The two components of @ on 0S, which approximate the 
components of Tu on that part of the boundary, are shown in Fig. 5.203 in Cartesian 
coordinates. Their parametric form for 0 < ¢ < 1, and the primary and secondary 
knots, can be seen in Fig. 5.204. These graphs indicate that the approximate density 
@ on OS; agrees very well with the exact one @ = Tu displayed in Figs. 5.199 and 
5.200. 
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Fig. 5.204 The components of 6[x(t]], 0 <+< 1 (parametric form). 


The components of @ on 0S2, which approximate those of u there, are graphed 
in Fig. 5.205. (We recall that two coefficients of two basis functions for both com- 
ponents have been modified to prevent jump discontinuities.) Their parametric form 
for 1 <t < 2, and the primary and secondary knots, can be seen in Fig. 5.206. These 
graphs indicate that the approximate density @ on OS» agrees very well with the 
exact one @[x] = u[x] x € OS, displayed in Figs. 5.197 and 5.198. 
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Fig. 5.205 The components of @[x], x € 0S2 (Cartesian coordinates). 
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Fig. 5.206 The components of @[x(t]], 1 <+ <2 (parametric form). 


Approximate solution. We use @, “, and @ in (5.76) to compute an approxima- 
tion @ to the exact solution u in S. The graphs of the two components of # are shown 
in Fig. 5.207. 
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Fig. 5.207 The components of di[x], x € S. 
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Error analysis. To visualize the difference between the approximate and exact 
densities, in Fig. 5.208 we graphed @ — @ relative to the maximum absolute value 
of ~ on OS. The same graphs in parametric form in Fig. 5.209, where the intervals 
0<t<1and1<t <2 correspond to 0S; and 0S, respectively, show that the 
relative error in @ on the entire boundary with the exception of the corners is less 
than 1%. 
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Fig. 5.209 The components of the relative error in @[x[t]] (parametric form). 


5.20 Remarks. (i) The direct method for the Dirichlet, Neumann, and mixed prob- 
lems appears to be well conditioned as well as computationally accurate. In general, 
comparing this method with the classic indirect method illustrated in the preceding 


sections, we conclude that the former has been, from a computational viewpoint, 
consistently superior to the latter. 


(ii) The set of coefficients {c1 jj, c2;,;} in (5.77), computed with the B-spline 
basis shown in Fig. 5.202, is 


{{{15.1, 17.4, 22.9, 34.3, 51., 75.1, 11L., 145., 168.}, 


£1.85, 5.02, 11.7, 14.6, 10.6, 5.31, 1.92, 0.877, 0.476}}, 
{{6.41, 6.64, 8.86, 11.8, 16.2, 21.3, 28.9, 33.9, 39.5}, 


{—95.3, 96.8, —91.6, 66.9, —41.5, —25., -16.8, —15., —14.8}}}. 
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